# Generate stress in 2D subduction **URL:** \ **Category:** PyLith\ **Tags:** generate-stress\ **Created:** [January 6, 2020, 4:47pm UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993 "2020-01-06T16:47:16Z")\ **Posts on this page:** 6\ **Page:** 1
**Author:** ![Lauriane](https://avatars.discourse-cdn.com/v4/letter/l/aca169/32.png) [@Lauriane](https://community.geodynamics.org/u/Lauriane)\ **Post date:** [January 6, 2020, 4:47pm UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993/1 "2020-01-06T16:47:16Z")
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 ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [January 6, 2020, 8:32pm UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993/2 "2020-01-06T20:32:40Z")
See the comment on line 19 of `generate_statedb.py`: ```auto 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: ```auto # 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) ``` ---
**Author:** ![Lauriane](https://avatars.discourse-cdn.com/v4/letter/l/aca169/32.png) [@Lauriane](https://community.geodynamics.org/u/Lauriane)\ **Post date:** [January 10, 2020, 1:13pm UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993/3 "2020-01-10T13:13:03Z")
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é ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [January 10, 2020, 8:08pm UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993/4 "2020-01-10T20:08:58Z")
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. ---
**Author:** ![Lauriane](https://avatars.discourse-cdn.com/v4/letter/l/aca169/32.png) [@Lauriane](https://community.geodynamics.org/u/Lauriane)\ **Post date:** [January 31, 2020, 10:15am UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993/5 "2020-01-31T10:15:32Z")
Dear Brad, Thank you for your reply. Is there already one script to generate the stress state with triangles? If yes, is it possible to share it with me? Thank you in advance, Best regards, Lauriane Baylé Le ven. 10 janv. 2020 à 21:19, Brad Aagaard via Computational Infrastructure for Geodynamics [geodynamics@discoursemail.com](mailto:geodynamics@discoursemail.com) a écrit : ---
**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\ **Post date:** [February 1, 2020, 11:22pm UTC](https://community.geodynamics.org/t/generate-stress-in-2d-subduction/993/6 "2020-02-01T23:22:02Z")
I don’t believe we have a script to generate the stress state information for triangles.