Skip to content

Commit 8c917ff

Browse files
dhyan1272dhyan1272
committed
StressDivergence added
1 parent 5ff1728 commit 8c917ff

8 files changed

Lines changed: 165 additions & 43 deletions

‎src/pmpo_MPMesh.cpp‎

Lines changed: 124 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1);
1111
void MPMesh::calculateStrain(){
1212
auto MPsPosition = p_MPs->getPositions();
1313
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
14+
auto MPsBasisGrads = p_MPs->getData<MPF_Basis_Grad_Vals>();
1415
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
1516
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
1617
//Mesh Fields
@@ -19,9 +20,6 @@ void MPMesh::calculateStrain(){
1920
auto gnomProjElmCenter = p_mesh->getMeshField<MeshF_ElmCenterGnomProj>();
2021
auto elm2VtxConn = p_mesh->getElm2VtxConn();
2122
auto velField = p_mesh->getMeshField<MeshF_Vel>();
22-
double radius = 1.0;
23-
if(p_mesh->getGeomType() == geom_spherical_surf)
24-
radius=p_mesh->getSphereRadius();
2523

2624
bool isRotated = p_mesh->getRotatedFlag();
2725

@@ -39,30 +37,28 @@ void MPMesh::calculateStrain(){
3937
computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
4038
auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL);
4139

42-
double basisByArea[maxVtxsPerElm] = {0.0};
43-
initArray(basisByArea,maxVtxsPerElm,0.0);
44-
double gradBasisByArea[2*maxVtxsPerElm] = {0.0};
45-
initArray(gradBasisByArea,maxVtxsPerElm,0.0);
46-
47-
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, radius, basisByArea, gradBasisByArea);
48-
4940
double v11 = 0.0;
5041
double v12 = 0.0;
5142
double v21 = 0.0;
5243
double v22 = 0.0;
5344
double uTanOverR = 0.0;
5445
double vTanOverR = 0.0;
46+
5547
for (int i = 0; i < numVtx; i++){
5648
int iVertex = elm2VtxConn(elm, i+1)-1;
57-
v11 = v11 + gradBasisByArea[i*2 + 0] * velField(iVertex, 0);
58-
v12 = v12 + gradBasisByArea[i*2 + 1] * velField(iVertex, 0);
59-
v21 = v21 + gradBasisByArea[i*2 + 0] * velField(iVertex, 1);
60-
v22 = v22 + gradBasisByArea[i*2 + 1] * velField(iVertex, 1);
61-
uTanOverR = uTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 0);
62-
vTanOverR = vTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1);
49+
v11 = v11 + MPsBasisGrads(mp, i*2 + 0) * velField(iVertex, 0);
50+
v12 = v12 + MPsBasisGrads(mp, i*2 + 1) * velField(iVertex, 0);
51+
v21 = v21 + MPsBasisGrads(mp, i*2 + 0) * velField(iVertex, 1);
52+
v22 = v22 + MPsBasisGrads(mp, i*2 + 1) * velField(iVertex, 1);
53+
uTanOverR = uTanOverR + MPsBasis(mp, i) * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 0);
54+
vTanOverR = vTanOverR + MPsBasis(mp, i) * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1);
6355
MPsStrainRate(mp, 0) = v11 - vTanOverR;
6456
MPsStrainRate(mp, 1) = v22;
6557
MPsStrainRate(mp, 2) = 0.5*(v12 + v21 + uTanOverR);
58+
//Debugging
59+
if(MPsAppID(mp)==0){
60+
printf("Strain Calc: iVertex %d vel field %.15e %.15e \n", iVertex, velField(iVertex, 0), velField(iVertex, 1));
61+
}
6662
}
6763
MPsStrainRate(mp, 0) = v11 - vTanOverR;
6864
MPsStrainRate(mp, 1) = v22;
@@ -79,32 +75,133 @@ void MPMesh::calculateStress(){
7975
auto dynamicTimeStep = p_mesh->getDynamicTimeStep();
8076
auto dampingTimescale = polyMPO::dampingTimescaleParameter * dynamicTimeStep;
8177
//MPFields
82-
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
83-
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
84-
auto MPsArea = p_MPs->getData<polyMPO::MPF_Area>();
78+
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
79+
auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
80+
auto MPsStress = p_MPs->getData<MPF_Stress>();
81+
auto MPsArea = p_MPs->getData<polyMPO::MPF_Area>();
8582
auto MPsIcePressure = p_MPs->getData<polyMPO::MPF_IcePressure>();
8683
auto MPsRepPressure = p_MPs->getData<polyMPO::MPF_ReplacementPressure>();
8784

8885
auto setMPStress = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
8986
if(mask){
87+
88+
Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
89+
Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 1));
90+
constitutive_evp(strain_rate, stress, MPsIcePressure(mp, 0), MPsRepPressure(mp, 0), MPsArea(mp, 0), elasticTimeStep, dampingTimescale);
91+
for (int m=0 ; m<3; m++)
92+
MPsStress(mp, m) = stress[m];
9093
//Debugging
9194
if(MPsAppID(mp)==0){
92-
printf("Strain %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
93-
printf("Mesh:%.15e %d, MP:%.15e %.15e %.15e\n",elasticTimeStep,solveStress(elm),MPsArea(mp,0),MPsIcePressure(mp,0),MPsRepPressure(mp,0));
94-
Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
95-
Vec3d stress(0, 0, 0);
96-
constitutive_evp(strain_rate, stress, MPsIcePressure(mp, 0), MPsRepPressure(mp, 0), MPsArea(mp, 0), elasticTimeStep, dampingTimescale);
95+
printf("Stress in GPU: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
9796
}
9897
}
9998
};
10099
p_MPs->parallel_for(setMPStress, "setMPStress");
101100
}
102101

102+
void MPMesh::calculateStressDivergence(){
103+
104+
//Mesh Information
105+
auto elm2VtxConn = p_mesh->getElm2VtxConn();
106+
int numVtx = p_mesh->getNumVertices();
107+
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
108+
int numVertices = p_mesh->getNumVertices();
109+
auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
110+
111+
//Material Points
112+
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
113+
auto weight = p_MPs->getData<MPF_Basis_Vals>();
114+
auto weight_grads = p_MPs->getData<MPF_Basis_Grad_Vals>();
115+
auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();
116+
auto MPsStress = p_MPs->getData<MPF_Stress>();
117+
auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
118+
119+
auto VtxCoeffs_new = this->precomputedVtxCoeffs_new;
120+
121+
//Earth Radius
122+
double radius = 1.0;
123+
if(p_mesh->getGeomType() == geom_spherical_surf)
124+
radius=p_mesh->getSphereRadius();
125+
126+
//Reconstructed the stress
127+
Kokkos::View<double*[3]> stress_rec("stress_rec", p_mesh->getNumVertices());
128+
Kokkos::View<double*> stress_divU("stress_divu", p_mesh->getNumVertices());
129+
Kokkos::View<double*> stress_divV("stress_divv", p_mesh->getNumVertices());
130+
Kokkos::View<double*> dSdX("dSdX", p_mesh->getNumVertices());
131+
Kokkos::View<double*> dSdY("dSdY", p_mesh->getNumVertices());
132+
133+
//Assemble fields for Stress Divergence
134+
auto stress_div = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
135+
if(mask) { //if material point is 'active'/'enabled'
136+
int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
137+
for(int i=0; i<nVtxE; i++){
138+
int vID = elm2VtxConn(elm,i+1)-1;
139+
double w_vtx=weight(mp,i);
140+
double CoordDiffs[vec4d_nEntries] = {1, (-vtxCoords(vID,0) + mpPositions(mp,0))/radius,
141+
(-vtxCoords(vID,1) + mpPositions(mp,1))/radius,
142+
(-vtxCoords(vID,2) + mpPositions(mp,2))/radius};
143+
144+
auto factor = w_vtx*(VtxCoeffs_new(vID,0, 0) + VtxCoeffs_new(vID,0, 1)*CoordDiffs[1] +
145+
VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] +
146+
VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]);
147+
148+
auto factor1 = (w_vtx/radius)*(VtxCoeffs_new(vID,1, 0) + VtxCoeffs_new(vID,1, 1)*CoordDiffs[1] +
149+
VtxCoeffs_new(vID,1, 2)*CoordDiffs[2] +
150+
VtxCoeffs_new(vID,1, 3)*CoordDiffs[3]);
151+
auto factor2 = (w_vtx/radius)*(VtxCoeffs_new(vID,2, 0) + VtxCoeffs_new(vID,2, 1)*CoordDiffs[1] +
152+
VtxCoeffs_new(vID,2, 2)*CoordDiffs[2] +
153+
VtxCoeffs_new(vID,2, 3)*CoordDiffs[3]);
154+
155+
for (int k=0; k<3; k++){
156+
auto val = factor*MPsStress(mp,k);
157+
Kokkos::atomic_add(&stress_rec(vID,k), val);
158+
}
159+
Kokkos::atomic_add(&stress_divU(vID), factor1 * MPsStress(mp, 0) + factor2 * MPsStress(mp, 2) -
160+
2 * tanLatVertexRotatedOverRadius(vID, 0) * factor * MPsStress(mp, 2));
161+
Kokkos::atomic_add(&stress_divV(vID), factor2 * MPsStress(mp, 1) + factor1*MPsStress(mp, 2) +
162+
factor * tanLatVertexRotatedOverRadius(vID, 0) * (MPsStress(mp, 0)-MPsStress(mp, 1)));
163+
Kokkos::atomic_add(&dSdX(vID), -factor * weight_grads(mp, i*2 + 0));
164+
Kokkos::atomic_add(&dSdY(vID), -factor * weight_grads(mp, i*2 + 1));
165+
}
166+
}
167+
};
168+
p_MPs->parallel_for(stress_div, "assembly");
169+
170+
//TODO
171+
//COMMUNICATE THE VERTEX FIELDS
172+
173+
//TODO put as mesh field
174+
Kokkos::View<doubleSclr_t*> stressDivergence("stressDivergence", p_mesh->getNumVertices());
175+
auto areaVertex = p_mesh->getMeshField<MeshF_DualTriangleArea>();
176+
177+
Kokkos::parallel_for("calculate_divergence", numVtx, KOKKOS_LAMBDA(const int vtx){
178+
double threshold = 0.15 / Kokkos::sqrt(areaVertex(vtx, 0));
179+
double valX = Kokkos::max(Kokkos::abs(dSdX(vtx)) - threshold, 0.0);
180+
double dSdX_filtered = Kokkos::copysign(valX, dSdX(vtx));
181+
double valY = Kokkos::max(Kokkos::abs(dSdY(vtx)) - threshold, 0.0);
182+
double dSdY_filtered = Kokkos::copysign(valY, dSdY(vtx));
183+
184+
stressDivergence(vtx, 0) = stress_divU(vtx) + dSdX_filtered * stress_rec(vtx, 0) + dSdY_filtered * stress_rec(vtx, 2);
185+
stressDivergence(vtx, 1) = stress_divV(vtx) + dSdY_filtered * stress_rec(vtx, 1) + dSdX_filtered * stress_rec(vtx, 2);
186+
//Debugging
187+
if (vtx >= 10 && vtx <= 11) {
188+
printf("Vtx %d Area %.15e ds: %.15e %.15e \n", vtx, areaVertex(vtx, 0), dSdX_filtered, dSdY_filtered);
189+
printf("Vtx %d Divergence %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
190+
}
191+
});
192+
193+
//TODO
194+
//COMMUNICATE THE VERTEX FIELDS
195+
196+
}
197+
198+
103199
void MPMesh::calcBasis() {
104200
assert(p_mesh->getGeomType() == geom_spherical_surf);
105201

106202
auto MPsPosition = p_MPs->getPositions();
107203
auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
204+
auto MPsBasisGrads = p_MPs->getData<MPF_Basis_Grad_Vals>();
108205
auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
109206

110207
auto elm2VtxConn = p_mesh->getElm2VtxConn();
@@ -141,7 +238,9 @@ void MPMesh::calcBasis() {
141238
wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, radius, basisByArea, gradBasisByArea);
142239

143240
for(int i=0; i<= numVtx; i++){
144-
MPsBasis(mp,i) = basisByArea[i];
241+
MPsBasis(mp, i) = basisByArea[i];
242+
MPsBasisGrads(mp, i*2+0) = gradBasisByArea[i*2 + 0];
243+
MPsBasisGrads(mp, i*2+1) = gradBasisByArea[i*2 + 1];
145244
}
146245

147246
//Old method where basis functions calculated using 3D Area

‎src/pmpo_MPMesh.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@ class MPMesh{
8888
void printVTP_mesh(int printVTPIndex);
8989
void calculateStrain();
9090
void calculateStress();
91+
void calculateStressDivergence();
9192
};
9293

9394
}//namespace polyMPO end

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -270,13 +270,23 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
270270
iBlockC[0] = blockC[0] * invM2D[0] + blockC[1] * invM2D[1] + blockC[2] * invM2D[2];
271271
iBlockC[1] = blockC[0] * invM2D[1] + blockC[1] * invM2D[3] + blockC[2] * invM2D[4];
272272
iBlockC[2] = blockC[0] * invM2D[2] + blockC[1] * invM2D[4] + blockC[2] * invM2D[5];
273-
iBlockC = iBlockC*invM11;
273+
iBlockC = iBlockC*invM11;
274274

275275
VtxCoeffs(vtx, 0, 0) = invM11 + invM11*iBlockC.dot(blockC);
276276
auto temp = -rotateScaleM.rightMultiply(iBlockC);
277277
VtxCoeffs(vtx, 0, 1) = temp[0];
278278
VtxCoeffs(vtx, 0, 2) = temp[1];
279279
VtxCoeffs(vtx, 0, 3) = temp[2];
280+
281+
VtxCoeffs(vtx, 1, 0) = -iBlockC[0];
282+
VtxCoeffs(vtx, 1, 1) = invM2D[0] * rotateScaleM(0, 0) + invM2D[1] * rotateScaleM(0, 1) + invM2D[2] * rotateScaleM(0, 2);
283+
VtxCoeffs(vtx, 1, 2) = invM2D[0] * rotateScaleM(1, 0) + invM2D[1] * rotateScaleM(1, 1) + invM2D[2] * rotateScaleM(1, 2);
284+
VtxCoeffs(vtx, 1, 3) = invM2D[0] * rotateScaleM(2, 0) + invM2D[1] * rotateScaleM(2, 1) + invM2D[2] * rotateScaleM(2, 2);
285+
286+
VtxCoeffs(vtx, 2, 0) = -iBlockC[1];
287+
VtxCoeffs(vtx, 2, 1) = invM2D[1] * rotateScaleM(0, 0) + invM2D[3] * rotateScaleM(0, 1) + invM2D[4] * rotateScaleM(0, 2);
288+
VtxCoeffs(vtx, 2, 2) = invM2D[1] * rotateScaleM(1, 0) + invM2D[3] * rotateScaleM(1, 1) + invM2D[4] * rotateScaleM(1, 2);
289+
VtxCoeffs(vtx, 2, 3) = invM2D[1] * rotateScaleM(2, 0) + invM2D[3] * rotateScaleM(2, 1) + invM2D[4] * rotateScaleM(2, 2);
280290
});
281291
this->precomputedVtxCoeffs_new = VtxCoeffs;
282292
}

‎src/pmpo_c.cpp‎

Lines changed: 18 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -698,7 +698,7 @@ void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int
698698
pumipic::RecordTime("PolyMPO_setIcePressure" + std::to_string(self), timer.seconds());
699699
}
700700

701-
void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){
701+
void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){
702702
Kokkos::Timer timer;
703703
checkMPMeshValid(p_mpmesh);
704704

@@ -713,18 +713,20 @@ void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, c
713713
//MP Data
714714
auto mpReplacementPressure = p_MPs->getData<polyMPO::MPF_ReplacementPressure>();
715715
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
716-
//Copy to device
717-
kkViewHostU<const double**> mpReplacementPressure_h(replacementPressureMPHost, nComps, numMPs);
718-
Kokkos::View<double**> mpReplacementPressure_d("mpIcePressureDevice", nComps, numMPs);
719-
Kokkos::deep_copy(mpReplacementPressure_d, mpReplacementPressure_h);
720-
//Set in PS
721-
auto setMPReplacementPressure = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
716+
//Copy to host
717+
Kokkos::View<double**> mpReplacementPressure_d("mpIcePressureDevice", nComps, numMPs);
718+
//From PS to array
719+
auto getMPReplacementPressure = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
722720
if(mask){
723-
mpReplacementPressure(mp,0) = mpReplacementPressure_d(0, mpAppID(mp));
721+
mpReplacementPressure_d(0, mpAppID(mp)) = mpReplacementPressure(mp,0);
724722
}
725723
};
726-
p_MPs->parallel_for(setMPReplacementPressure, "setReplacementPressure");
727-
pumipic::RecordTime("PolyMPO_setReplacementPressure" + std::to_string(self), timer.seconds());
724+
p_MPs->parallel_for(getMPReplacementPressure, "setReplacementPressure");
725+
//GPU to CPU
726+
kkDbl2dViewHostU arrayHost(replacementPressureMPHost, nComps, numMPs);
727+
Kokkos::deep_copy(arrayHost, mpReplacementPressure_d);
728+
pumipic::RecordTime("PolyMPO_getReplacementPressure" + std::to_string(self), timer.seconds());
729+
728730
}
729731

730732
void polympo_startMeshFill_f(MPMesh_ptr p_mpmesh){
@@ -1327,6 +1329,12 @@ void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* ar
13271329
Kokkos::deep_copy(solveStress, h_solveStress);
13281330
}
13291331

1332+
void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh){
1333+
//chech validity
1334+
checkMPMeshValid(p_mpmesh);
1335+
((polyMPO::MPMesh*)p_mpmesh) -> calculateStressDivergence();
1336+
}
1337+
13301338
//Advection Calcualtions
13311339
void polympo_push_f(MPMesh_ptr p_mpmesh){
13321340
checkMPMeshValid(p_mpmesh);

‎src/pmpo_c.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@ void polympo_setMPStress_f(MPMesh_ptr p_mpmesh);
5252
void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost);
5353
void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost);
5454
void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost);
55-
void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost);
55+
void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost);
5656

5757
//Mesh info
5858
void polympo_startMeshFill_f(MPMesh_ptr p_mpmesh);
@@ -101,6 +101,7 @@ void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int n
101101
void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeStep);
102102
void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep);
103103
void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array);
104+
void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh);
104105

105106
// Advection calculations
106107
void polympo_push_f(MPMesh_ptr p_mpmesh);

‎src/pmpo_const_relation.hpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,6 @@ static constexpr double dampingTimescaleParameter = 0.36;
1414
KOKKOS_INLINE_FUNCTION
1515
void constitutive_evp(const Vec3d& strain, Vec3d& stress, const double& icePressure, double& replacementPressure,
1616
const double& areaMP, const double& dtElastic, const double& dampingTimescale){
17-
printf("In the Constitutive model\n");
1817
auto strainDivergence = strain[0] + strain[1];
1918
auto strainTension = strain[0] - strain[1];
2019
auto strainShearing = 2*strain[2];
@@ -36,8 +35,6 @@ void constitutive_evp(const Vec3d& strain, Vec3d& stress, const double& icePress
3635

3736
stress[0] = 0.5 * (stress1 + stress2);
3837
stress[1] = 0.5 * (stress1 - stress2);
39-
40-
printf("Stress %.15e %.15e %.15e \n", stress[0], stress[1], stress[2]);
4138
}
4239

4340
}

‎src/pmpo_fortran.f90‎

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -393,8 +393,8 @@ subroutine polympo_setIcePressureMP(mpMesh, nComps, numMPs, array) &
393393
type(c_ptr), value :: array
394394
end subroutine
395395

396-
subroutine polympo_setReplacementPressureMP(mpMesh, nComps, numMPs, array) &
397-
bind(C, NAME='polympo_setReplacementPressureMP_f')
396+
subroutine polympo_getReplacementPressureMP(mpMesh, nComps, numMPs, array) &
397+
bind(C, NAME='polympo_getReplacementPressureMP_f')
398398
use :: iso_c_binding
399399
type(c_ptr), value :: mpMesh
400400
integer(c_int), value :: nComps, numMPs
@@ -950,6 +950,12 @@ subroutine polympo_setSolveStressMesh(mpMesh, nCells, array) &
950950
type(c_ptr), value :: array
951951
end subroutine
952952

953+
subroutine polympo_calculateStressDivergence(mpMesh) &
954+
bind(C, NAME='polympo_calculateStressDivergence_f')
955+
use :: iso_c_binding
956+
type(c_ptr), value :: mpMesh
957+
end subroutine
958+
953959
!---------------------------------------------------------------------------
954960
!> @brief calculate the MPs from given mesh vertices rotational latitude
955961
!> longitude, update the MP slices

‎src/pmpo_materialPoints.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@ template <> struct mpSliceToMeshField < MPF_Vel > { using type =
6868
template <> struct mpSliceToMeshField < MPF_Rot_Lat_Lon_Incr > { using type = vec2d_t; };
6969
template <> struct mpSliceToMeshField < MPF_Vel_Incr > { using type = vec2d_t; };
7070
template <> struct mpSliceToMeshField < MPF_Strain_Rate > { using type = double[3]; };
71-
template <> struct mpSliceToMeshField < MPF_Stress > { using type = double[6]; };
71+
template <> struct mpSliceToMeshField < MPF_Stress > { using type = double[3]; };
7272
template <> struct mpSliceToMeshField < MPF_Stress_Div > { using type = vec3d_t; };
7373
template <> struct mpSliceToMeshField < MPF_Shear_Traction > { using type = vec3d_t; };
7474
template <> struct mpSliceToMeshField < MPF_Constv_Mdl_Param > { using type = double[12]; };

0 commit comments

Comments
 (0)