 # Generate stress in 2D subduction

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
File “./generate_statedb.py”, line 42, in calcQuadCoords
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?

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

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)
``````