<div dir="ltr">Thanks Dave and Mathew. I will keep looking for the bug and here is the full piece of code.<div><br></div><div>PetscInt kel_tet10CX(PetscScalar *X, PetscScalar *Y, PetscScalar *Z, PetscScalar *ke, PetscInt elmat, PetscScalar dmg, PetscInt el_has_orient, PetscScalar rot[3][3], PetscScalar temperature, PetscScalar *thermforce, PetscInt damflag, PetscScalar Cmat[6][6], PetscScalar alphavec[6]){<br> /////////PETSc declarations <br> PetscErrorCode ierr;<br> PetscInt ii,jj,kk,k1,ll,l1,i,j,k,l,i1,j1,jj1,ii1,iflag;<br> PetscInt nope=10; //// nodes on the element<br> PetscScalar alpha2[4][4],alpha2unrot[4][4],C[6][6],C4[4][4][4][4],C4unrot[4][4][4][4],eta,etas[9],GP[4],s[100][100],shp[10][4],shpj[4][11],w[4][4],weight,wts[4],xi,xis[4],xl[10][3],xsj,xsjj,zeta,zetas[4];<br><br> /////////////Prototypes<br> PetscInt constit(PetscInt elmat,PetscScalar C4[4][4][4][4], PetscScalar alpha2[4][4], PetscScalar dmg, PetscScalar temperature, PetscInt damflag, PetscScalar Cmat[6][6], PetscScalar alphavec[6]);<br> void shape10tet_(PetscScalar *xi, PetscScalar *eta, PetscScalar *zeta, PetscScalar xl[4][11], PetscScalar *xsj, PetscScalar shp[5][11], PetscInt *iflag); // you have to dimension the array arguments to the size Fortran sees them<br><br> // printf("\n Entered SPAF element Tet10 \n");<br> /////////////Code<br> if (el_has_orient == 1){ierr = constit(elmat,C4unrot,alpha2unrot,dmg,temperature,damflag,Cmat,alphavec);ierr = rot4(rot,C4unrot,C4);ierr = rot2(rot,alpha2unrot,alpha2);}<br> else {ierr = constit(elmat,C4,alpha2,dmg,temperature,damflag,Cmat,alphavec);}<br><br> // printf("\n Came back from spaf materials to SPAF element Tet10 \n");<br> // ierr = PetscPrintf(PETSC_COMM_SELF,"In tet10 \n");CHKERRQ(ierr);<br><br><br> // Gauss points and weights<br> xis[0] = 0.138196601125011;<br> xis[1] = 0.585410196624968;<br> xis[2] = 0.138196601125011;<br> xis[3] = 0.138196601125011;<br> etas[0] = 0.138196601125011;<br> etas[1] = 0.138196601125011;<br> etas[2] = 0.585410196624968;<br> etas[3] = 0.138196601125011;<br> zetas[0] = 0.138196601125011;<br> zetas[1] = 0.138196601125011;<br> zetas[2] = 0.138196601125011;<br> zetas[3] = 0.585410196624968;<br> wts[0] = 0.04166666666667;<br> wts[1] = 0.04166666666667;<br> wts[2] = 0.04166666666667;<br> wts[3] = 0.04166666666667;<br><br> // assemble nodal coordinates into xl<br> // when fortran gets it, it will be transposed and the indices increased by 1<br> for (ii=0; ii<nope; ii++)<br> {<br> xl[ii][0] = X[ii];//fprintf(outfile,"X[%i] = %f\n",ii,X[ii]);<br> xl[ii][1] = Y[ii];//fprintf(outfile,"Y[%i] = %f\n",ii,Y[ii]);<br> xl[ii][2] = Z[ii];//fprintf(outfile,"Z[%i] = %f\n",ii,Z[ii]);<br> }<br><br> // ke<br> memset(ke, 0, sizeof(ke[0])*nope*3*nope*3); // why doesn't compiler know the real size of ke? <br> memset(s, 0, sizeof(s));<br> iflag = 3; /// tells shape routine to calculate shape functions, derivatives and det J<br><br> for (kk=0; kk<4; kk++) // Loop 3<br> {<br> weight = wts[kk];<br> // Integration point<br> xi = xis[kk]; <br> eta = etas[kk]; <br> zeta = zetas[kk];<br> shape10tet_(&xi,&eta,&zeta,xl,&xsj,shp,&iflag);<br> if (xsj<0){fprintf(outfile,"Negative det J on TET10 element!! (Integration point %i)\n",kk);}<br> xsjj = sqrt(xsj);<br> for (i1=1; i1<nope+1; i1++) // Loop 4a<br> {<br> shpj[1][i1] = shp[i1-1][0]*xsjj; // shape function derivative wrt xi<br> shpj[2][i1] = shp[i1-1][1]*xsjj; // shape function derivative wrt eta<br> shpj[3][i1] = shp[i1-1][2]*xsjj; // shape function derivative wrt zeta<br> shpj[4][i1] = shp[i1-1][3]*xsj; // shape function itself<br> // printf("\n shp[1][il] %f %f %f %f \n",shpj[1][i1],shpj[2][i1],shpj[3][i1],shpj[4][i1]);<br> } // Loop 4a<br> jj1=1;<br> for (j=1; j<nope+1; j++) // Loop 4b<br> {<br> ii1=1;<br> for (i=1; i<nope+1; i++) // Loop 5<br> {<br> for (i1=1; i1<4; i1++) // Loop 6a<br> {<br> for (j1=1; j1<4; j1++) // Loop 7a<br> {<br> w[i1][j1] = shpj[i1][i]*shpj[j1][j];<br><br> } // Loop 7a<br> } // Loop 6a<br> for (i1=1; i1<4; i1++) // Loop 6b<br> {<br> for (j1=1; j1<4; j1++) // Loop 7b<br> {<br> for (k1=1; k1<4; k1++) // Loop 8<br> {<br> for (l1=1; l1<4; l1++) // Loop 9<br> {<br> s[ii1+i1-1][jj1+j1-1] = s[ii1+i1-1][jj1+j1-1]+C4[i1][k1][j1][l1]*w[k1][l1]*weight;<br> ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"C4 %f \n",C4[i1][k1][j1][l1]);CHKERRQ(ierr);<br> ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"w %f \n",w[k1][l1]);CHKERRQ(ierr);<br> ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"weight %f \n",weight);CHKERRQ(ierr);<br> } // Loop 9<br> } // Loop 8<br> } // Loop 7b<br> } // Loop 6b<br> ii1 = ii1 + 3;<br> } // Loop 5<br> jj1 = jj1 + 3;<br> } // Loop 4b<br> } // Loop 3<br> // ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"C4 %f \n",C4[1][1][1][1]);CHKERRQ(ierr);<br><br><br> for (ii=1; ii<nope*3+1; ii++)<br> {<br> for (jj=1; jj<nope*3+1; jj++)<br> { ke[(ii-1)*nope*3+(jj-1)] = s[ii][jj];<br> // printf("\n ke[(%d-1)*nope*3+(%d-1)] = %f \n", ii,jj, ke[(ii-1)*nope*3+(jj-1)] );<br> }<br> }<br><br> // ierr = PetscPrintf(PETSC_COMM_SELF,"exiting spaf elemnts tet10 %f\n",ke[3200]);CHKERRQ(ierr);<br><br> return 0;<br><br>}<br></div><div><br></div><div>Thanks in advance</div><div>Kaushik</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Apr 21, 2020 at 10:51 AM Dave May <<a href="mailto:dave.mayhem23@gmail.com">dave.mayhem23@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><br></div><div><br><div class="gmail_quote"><div dir="ltr">On Tue 21. Apr 2020 at 16:47, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">You are overwriting memory somewhere. The prints just move it around. I suggest running with valgrind.<div></div></div></blockquote><div dir="auto"><br></div><div dir="auto">Matt is right. However, judging by the code snippet I bet all the arrays in question are statically allocated, thus valgrind may be of somewhat limited use.</div><div dir="auto"><br></div><div dir="auto">If you send the entire function, or all the related pieces of code, someone in the list might spot the error.</div><div dir="auto"><br></div><div dir="auto">Thanks</div><div dir="auto">Dave</div><div dir="auto"><br></div><div dir="auto"><br></div><div dir="auto"><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Apr 21, 2020 at 10:44 AM Kaushik Vijaykumar <<a href="mailto:kaushikv318@gmail.com" target="_blank">kaushikv318@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Hello group,<div><br></div><div>I have been trying to navigate a weird error that I have found in my FEA code that I am developing using PetSc. The error occurs when, i execute a call to stiffness generation of a Tetrahederal element. The code returns an memory error if I don't print the following statements in the function, (see the ierr print statements below):</div><div><br></div><div>for (i1=1; i1<4; i1++) // Loop 6b<br> {<br> for (j1=1; j1<4; j1++) // Loop 7b<br> {<br> for (k1=1; k1<4; k1++) // Loop 8<br> {<br> for (l1=1; l1<4; l1++) // Loop 9<br> {<br> s[ii1+i1-1][jj1+j1-1] = s[ii1+i1-1][jj1+j1-1]+C4[i1][k1][j1][l1]*w[k1][l1]*weight;<br> ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"C4 %f \n",C4[i1][k1][j1][l1]);CHKERRQ(ierr);<br> ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"w %f \n",w[k1][l1]);CHKERRQ(ierr);<br> ierr = PetscFPrintf(PETSC_COMM_WORLD,outfile,"weight %f \n",weight);CHKERRQ(ierr);<br> } // Loop 9<br> } // Loop 8<br> } // Loop 7b<br> } // Loop 6b<br></div><div><br></div><div><br></div><div>Any help on this is really appreciated.</div><div><br></div><div>Thanks</div><div>Kaushik</div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div>
</blockquote></div></div>
</blockquote></div>