Generate stress in 2D subduction

#1

Hi,

I am trying to initial stress in a 2D subduction zone with a triangle mesh. I put the gravity (using the example) in my subduction (same example than in the manual) and now I am trying to do the step8. I already re-run the step3 with commenting out the CellFilterAvg settings and run pylith gravity_initstress.cfg gravity_vardensity.cfg. But when I did ./generate_statedb.py, I have this error:
“Traceback (most recent call last):
File “./generate_statedb.py”, line 63, in
quadCoords = calcQuadCoords(vertices, cells, qpts)
File “./generate_statedb.py”, line 42, in calcQuadCoords
quadCoords[:,:,iDim] = numpy.dot(cellCoords[:,:,iDim], qpts.transpose())
ValueError: shapes (523,3) and (4,4) not aligned: 3 (dim 1) != 4 (dim 0)”

I understood that in the generate_statedb.py file the qpts array are the basis function evaluated at the quadrature points but I still don’t understand how to calculate this value? Furthermore, the example of 2D gravity is for lines and my model is for triangle so I need to say in generate_statedb.py the good size. But I don’t understand how?
Can you help me?

Thank you in advance,

Lauriane Baylé
Phd Student Institut des Sciences de la Terre de Paris
Sorbonne University

#2

See the comment on line 19 of generate_statedb.py:

 19 # Basis functions for quad4 cell evaluated at quadrature points. Use
 20 # to compute coordinate of quadrature points in each cell from
 21 # coordinates of vertices. Note the order must correspond to the order
 22 # of the data at the quadrature points in the output.
 23 qpts = numpy.array([[ 0.62200847,  0.16666667,  0.0446582,   0.16666667],
 24                     [ 0.16666667,  0.62200847,  0.16666667,  0.0446582 ],
 25                     [ 0.16666667,  0.0446582,   0.16666667,  0.62200847],
 26                     [ 0.0446582,   0.16666667,  0.62200847,  0.16666667]], dtype=numpy.float64)

For a mesh with triangle cells, you will need to update this qpts array (for quad cells) with the coordinates of the quadrature points in the reference triangular cell. Assuming you are using one-point quadrature, you will have:

# Basis functions for tri3 cell evaluated at quadrature points. Use
# to compute coordinate of quadrature points in each cell from
# coordinates of vertices. Note the order must correspond to the order
# of the data at the quadrature points in the output.
qpts = numpy.array([[0.33333333, 0.33333333, 0.3333333]], dtype=numpy.float64)
#3

Thank you for your reply.
I am wondering if the qpts are the barycenter of the triangle? Furthermore, is it possible to have a calcTriCoords and so TriCoords?

Thank you,
Best regards,

Lauriane Baylé

#4

For one-point quadrature the point is in the center of the triangle (note that the basis functions are all equal at the quadrature point).

Yes, you can add functions to the script as necessary.