"Ferrand, Jesus A." <ferra...@my.erau.edu> writes: > Dear PETSc Team: > > I have a question about DM and PetscSection. Say I import a mesh (for FEM > purposes) and create a DMPlex for it. I then use PetscSections to set degrees > of freedom per "point" (by point I mean vertices, lines, faces, and cells). I > then use PetscSectionGetStorageSize() to get the size of the global stiffness > matrix (K) needed for my FEM problem. One last detail, this K I populate > inside a rather large loop using an element stiffness matrix function of my > own. Instead of using DMCreateMatrix(), I manually created a Mat using > MatCreate(), MatSetType(), MatSetSizes(), and MatSetUp(). I come to find that > said loop is painfully slow when I use the manually created matrix, but 20x > faster when I use the Mat coming out of DMCreateMatrix().
The sparse matrix hasn't been preallocated, which forces the data structure to do a lot of copies (as bad as O(n^2) complexity). DMCreateMatrix() preallocates for you. https://petsc.org/release/docs/manual/performance/#memory-allocation-for-sparse-matrix-assembly https://petsc.org/release/docs/manual/mat/#sec-matsparse > My question is then: Is the manual Mat a noob mistake and is it somehow > creating a memory leak with K? Just in case it's something else I'm attaching > the code. The loop that populates K is between lines 221 and 278. Anything > related to DM, DMPlex, and PetscSection is between lines 117 and 180. > > Machine Type: HP Laptop > C-compiler: Gnu C > OS: Ubuntu 20.04 > PETSc version: 3.16.0 > MPI Implementation: MPICH > > Hope you all had a Merry Christmas and that you will have a happy and > productive New Year. :D > > > Sincerely: > > J.A. Ferrand > > Embry-Riddle Aeronautical University - Daytona Beach FL > > M.Sc. Aerospace Engineering | May 2022 > > B.Sc. Aerospace Engineering > > B.Sc. Computational Mathematics > > > > Sigma Gamma Tau > > Tau Beta Pi > > Honors Program > > > > Phone: (386)-843-1829 > > Email(s): ferra...@my.erau.edu > > jesus.ferr...@gmail.com > //REFERENCE: > https://github.com/FreeFem/FreeFem-sources/blob/master/plugin/mpi/PETSc-code.hpp > #include <petsc.h> > static char help[] = "Imports a Gmsh mesh with boundary conditions and solves > the elasticity equation.\n" > "Option prefix = opt_.\n"; > > struct preKE{//Preallocation before computing KE > Mat matB, > matBTCB; > //matKE; > PetscInt x_insert[3], > y_insert[3], > z_insert[3], > m,//Looping variables. > sizeKE,//size of the element stiffness matrix. > N,//Number of nodes in element. > x_in,y_in,z_in; //LI to index B matrix. > PetscReal J[3][3],//Jacobian matrix. > invJ[3][3],//Inverse of the Jacobian matrix. > detJ,//Determinant of the Jacobian. > dX[3], > dY[3], > dZ[3], > minor00, > minor01, > minor02,//Determinants of minors in a 3x3 matrix. > dPsidX, dPsidY, dPsidZ,//Shape function derivatives w.r.t global > coordinates. > weight,//Multiplier of quadrature weights. > *dPsidXi,//Derivatives of shape functions w.r.t. Xi. > *dPsidEta,//Derivatives of shape functions w.r.t. Eta. > *dPsidZeta;//Derivatives of shape functions w.r.t Zeta. > PetscErrorCode ierr; > };//end struct. > > //Function declarations. > extern PetscErrorCode tetra4(PetscScalar*, PetscScalar*, PetscScalar*,struct > preKE*, Mat*, Mat*); > extern PetscErrorCode ConstitutiveMatrix(Mat*,const char*,PetscInt); > extern PetscErrorCode InitializeKEpreallocation(struct preKE*,const char*); > > PetscErrorCode PetscViewerVTKWriteFunction(PetscObject vec,PetscViewer > viewer){ > PetscErrorCode ierr; > ierr = VecView((Vec)vec,viewer); CHKERRQ(ierr); > return ierr; > } > > > > > int main(int argc, char **args){ > //DEFINITIONS OF PETSC's DMPLEX LINGO: > //POINT: A topology element (cell, face, edge, or vertex). > //CHART: It an interval from 0 to the number of "points." (the range of > admissible linear indices) > //STRATUM: A subset of the "chart" which corresponds to all "points" at a > given "level." > //LEVEL: This is either a "depth" or a "height". > //HEIGHT: Dimensionality of an element measured from 0D to 3D. Heights: > cell = 0, face = 1, edge = 2, vertex = 3. > //DEPTH: Dimensionality of an element measured from 3D to 0D. Depths: cell > = 3, face = 2, edge = 1, vertex = 0; > //CLOSURE: *of an element is the collection of all other elements that > define it.I.e., the closure of a surface is the collection of vertices and > edges that make it up. > //STAR: > //STANDARD LABELS: These are default tags that DMPlex has for its topology. > ("depth") > PetscErrorCode ierr;//Error tracking variable. > DM dm;//Distributed memory object (useful for managing grids.) > DMLabel physicalgroups;//Identifies user-specified tags in gmsh (to impose > BC's). > DMPolytopeType celltype;//When looping through cells, determines its type > (tetrahedron, pyramid, hexahedron, etc.) > PetscSection s; > KSP ksp;//Krylov Sub-Space (linear solver object) > Mat K,//Global stiffness matrix (Square, assume unsymmetric). > KE,//Element stiffness matrix (Square, assume unsymmetric). > matC;//Constitutive matrix. > Vec XYZ,//Coordinate vector, contains spatial locations of mesh's vertices > (NOTE: This vector self-destroys!). > U,//Displacement vector. > F;//Load Vector. > PetscViewer XYZviewer,//Viewer object to output mesh to ASCII format. > XYZpUviewer; //Viewer object to output displacements to ASCII > format. > PetscBool interpolate = PETSC_TRUE,//Instructs Gmsh importer whether to > generate faces and edges (Needed when using P2 or higher elements). > useCone = PETSC_TRUE,//Instructs "DMPlexGetTransitiveClosure()" > whether to extract the closure or the star. > dirichletBC = PETSC_FALSE,//For use when assembling the K matrix. > neumannBC = PETSC_FALSE,//For use when assembling the F vector. > saveASCII = PETSC_FALSE,//Whether to save results in ASCII format. > saveVTK = PETSC_FALSE;//Whether to save results as VTK format. > PetscInt nc,//number of cells. (PETSc lingo for "elements") > nv,//number of vertices. (PETSc lingo for "nodes") > nf,//number of faces. (PETSc lingo for "surfaces") > ne,//number of edges. (PETSc lingo for "lines") > pStart,//starting LI of global elements. > pEnd,//ending LI of all elements. > cStart,//starting LI for cells global arrangement. > cEnd,//ending LI for cells in global arrangement. > vStart,//starting LI for vertices in global arrangement. > vEnd,//ending LI for vertices in global arrangement. > fStart,//starting LI for faces in global arrangement. > fEnd,//ending LI for faces in global arrangement. > eStart,//starting LI for edges in global arrangement. > eEnd,//ending LI for edges in global arrangement. > sizeK,//Size of the element stiffness matrix. > ii,jj,kk,//Dedicated looping variables. > indexXYZ,//Variable to access the elements of XYZ vector. > indexK,//Variable to access the elements of the U and F vectors > (can reference rows and colums of K matrix.) > *closure = PETSC_NULL,//Pointer to the closure elements of a cell. > size_closure,//Size of the closure of a cell. > dim,//Dimension of the mesh. > //*edof,//Linear indices of dof's inside the K matrix. > dof = 3,//Degrees of freedom per node. > cells=0, edges=0, vertices=0, faces=0,//Topology counters when > looping through cells. > pinXcode=10, pinZcode=11,forceZcode=12;//Gmsh codes to extract > relevant "Face Sets." > PetscReal //*x_el,//Pointer to a vector that will store the x-coordinates > of an element's vertices. > //*y_el,//Pointer to a vector that will store the y-coordinates > of an element's vertices. > //*z_el,//Pointer to a vector that will store the z-coordinates > of an element's vertices. > *xyz_el,//Pointer to xyz array in the XYZ vector. > traction = -10, > *KEdata, > t1,t2; //time keepers. > const char *gmshfile = "TopOptmeshfine2.msh";//Name of gmsh file to import. > > ierr = PetscInitialize(&argc,&args,NULL,help); if(ierr) return ierr; //And > the machine shall work.... > > //MESH > IMPORT================================================================= > //IMPORTANT NOTE: Gmsh only creates "cells" and "vertices," it does not > create the "faces" or "edges." > //Gmsh probably can generate them, must figure out how to. > t1 = MPI_Wtime(); > ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,gmshfile,interpolate,&dm); > CHKERRQ(ierr);//Read Gmsh file and generate the DMPlex. > ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr);//1-D, 2-D, or 3-D > ierr = DMPlexGetChart(dm, &pStart, &pEnd); CHKERRQ(ierr);//Extracts linear > indices of cells, vertices, faces, and edges. > ierr = DMGetCoordinatesLocal(dm,&XYZ); CHKERRQ(ierr);//Extracts coordinates > from mesh.(Contiguous storage: [x0,y0,z0,x1,y1,z1,...]) > ierr = VecGetArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's > data. > t2 = MPI_Wtime(); > PetscPrintf(PETSC_COMM_WORLD,"Mesh Import time: %10f\n",t2-t1); > DMView(dm,PETSC_VIEWER_STDOUT_WORLD); > > //IMPORTANT NOTE: PETSc assumes that vertex-cell meshes are 2D even if they > were 3D, so its ordering changes. > //Cells remain at height 0, but vertices move to height 1 from height 3. To > prevent this from becoming an issue > //the "interpolate" variable is set to PETSC_TRUE to tell the mesh importer > to generate faces and edges. > //PETSc, therefore, technically does additional meshing. Gotta figure out > how to get this from Gmsh directly. > ierr = DMPlexGetDepthStratum(dm,3, &cStart, &cEnd);//Get LI of cells. > ierr = DMPlexGetDepthStratum(dm,2, &fStart, &fEnd);//Get LI of faces > ierr = DMPlexGetDepthStratum(dm,1, &eStart, &eEnd);//Get LI of edges. > ierr = DMPlexGetDepthStratum(dm,0, &vStart, &vEnd);//Get LI of vertices. > ierr = DMGetStratumSize(dm,"depth", 3, &nc);//Get number of cells. > ierr = DMGetStratumSize(dm,"depth", 2, &nf);//Get number of faces. > ierr = DMGetStratumSize(dm,"depth", 1, &ne);//Get number of edges. > ierr = DMGetStratumSize(dm,"depth", 0, &nv);//Get number of vertices. > /* > PetscPrintf(PETSC_COMM_WORLD,"global start = %10d\t global end = > %10d\n",pStart,pEnd); > PetscPrintf(PETSC_COMM_WORLD,"#cells = %10d\t i = %10d\t i < > %10d\n",nc,cStart,cEnd); > PetscPrintf(PETSC_COMM_WORLD,"#faces = %10d\t i = %10d\t i < > %10d\n",nf,fStart,fEnd); > PetscPrintf(PETSC_COMM_WORLD,"#edges = %10d\t i = %10d\t i < > %10d\n",ne,eStart,eEnd); > PetscPrintf(PETSC_COMM_WORLD,"#vertices = %10d\t i = %10d\t i < > %10d\n",nv,vStart,vEnd); > */ > //MESH > IMPORT================================================================= > > //NOTE: This section extremely hardcoded right now. > //Current setup would only support P1 meshes. > //MEMORY ALLOCATION > ========================================================== > ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s); CHKERRQ(ierr); > //The chart is akin to a contiguous memory storage allocation. Each chart > entry is associated > //with a "thing," could be a vertex, face, cell, or edge, or anything else. > ierr = PetscSectionSetChart(s, pStart, pEnd); CHKERRQ(ierr); > //For each "thing" in the chart, additional room can be made. This is > helpful for associating > //nodes to multiple degrees of freedom. These commands help associate nodes > with > for(ii = cStart; ii < cEnd; ii++){//Cell loop. > ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no > dof's associated with cells. > for(ii = fStart; ii < fEnd; ii++){//Face loop. > ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no > dof's associated with faces. > for(ii = vStart; ii < vEnd; ii++){//Vertex loop. > ierr = PetscSectionSetDof(s, ii, dof);CHKERRQ(ierr);}//Sets x, y, and z > displacements as dofs. > for(ii = eStart; ii < eEnd; ii++){//Edge loop > ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no > dof's associated with edges. > ierr = PetscSectionSetUp(s); CHKERRQ(ierr); > ierr = PetscSectionGetStorageSize(s,&sizeK);CHKERRQ(ierr);//Determine the > size of the global stiffness matrix. > ierr = DMSetLocalSection(dm,s); CHKERRQ(ierr);//Associate the PetscSection > with the DM object. > //PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec) > //ierr = DMCreateGlobalVector(dm,&U); CHKERRQ(ierr); > PetscSectionDestroy(&s); > //PetscPrintf(PETSC_COMM_WORLD,"sizeK = %10d\n",sizeK); > > //OBJECT > SETUP================================================================ > //Global stiffness matrix. > //PetscErrorCode DMCreateMatrix(DM dm,Mat *mat) > > //This makes the loop fast. > ierr = DMCreateMatrix(dm,&K); > > //This makes the loop uber slow. > //ierr = MatCreate(PETSC_COMM_WORLD,&K); CHKERRQ(ierr); > //ierr = MatSetType(K,MATAIJ); CHKERRQ(ierr);// Global stiffness matrix set > to some sparse type. > //ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,sizeK,sizeK); > CHKERRQ(ierr); > //ierr = MatSetUp(K); CHKERRQ(ierr); > > //Displacement vector. > ierr = VecCreate(PETSC_COMM_WORLD,&U); CHKERRQ(ierr); > ierr = VecSetType(U,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix > set to some sparse type. > ierr = VecSetSizes(U,PETSC_DECIDE,sizeK); CHKERRQ(ierr); > > //Load vector. > ierr = VecCreate(PETSC_COMM_WORLD,&F); CHKERRQ(ierr); > ierr = VecSetType(F,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix > set to some sparse type. > ierr = VecSetSizes(F,PETSC_DECIDE,sizeK); CHKERRQ(ierr); > //OBJECT > SETUP================================================================ > > //WARNING: This loop is currently hardcoded for P1 elements only! Must > Figure > //out a clever way to modify to accomodate Pn (n>1) elements. > > //BEGIN GLOBAL STIFFNESS MATRIX > BUILDER======================================= > t1 = MPI_Wtime(); > > > //PREALLOCATIONS============================================================== > ierr = ConstitutiveMatrix(&matC,"isotropic",0); CHKERRQ(ierr); > struct preKE preKEtetra4; > ierr = InitializeKEpreallocation(&preKEtetra4,"tetra4"); CHKERRQ(ierr); > ierr = MatCreate(PETSC_COMM_WORLD,&KE); CHKERRQ(ierr); //SEQUENTIAL > ierr = MatSetSizes(KE,PETSC_DECIDE,PETSC_DECIDE,12,12); CHKERRQ(ierr); > //SEQUENTIAL > ierr = MatSetType(KE,MATDENSE); CHKERRQ(ierr); //SEQUENTIAL > ierr = MatSetUp(KE); CHKERRQ(ierr); > PetscReal x_tetra4[4], y_tetra4[4],z_tetra4[4], > x_hex8[8], y_hex8[8],z_hex8[8], > *x,*y,*z; > PetscInt *EDOF,edof_tetra4[12],edof_hex8[24]; > DMPolytopeType previous = DM_POLYTOPE_UNKNOWN; > > //PREALLOCATIONS============================================================== > > > > for(ii=cStart;ii<cEnd;ii++){//loop through cells. > ierr = DMPlexGetTransitiveClosure(dm, ii, useCone, &size_closure, > &closure); CHKERRQ(ierr); > ierr = DMPlexGetCellType(dm, ii, &celltype); CHKERRQ(ierr); > //IMPORTANT NOTE: MOST OF THIS LOOP SHOULD BE INCLUDED IN THE KE3D > function. > if(previous != celltype){ > //PetscPrintf(PETSC_COMM_WORLD,"run \n"); > if(celltype == DM_POLYTOPE_TETRAHEDRON){ > x = x_tetra4; > y = y_tetra4; > z = z_tetra4; > EDOF = edof_tetra4; > }//end if. > else if(celltype == DM_POLYTOPE_HEXAHEDRON){ > x = x_hex8; > y = y_hex8; > z = z_hex8; > EDOF = edof_hex8; > }//end else if. > } > previous = celltype; > > //PetscPrintf(PETSC_COMM_WORLD,"Cell # %4i\t",ii); > cells=0; > edges=0; > vertices=0; > faces=0; > kk = 0; > for(jj=0;jj<(2*size_closure);jj+=2){//Scan the closure of the current > cell. > //Use information from the DM's strata to determine composition of > cell_ii. > if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices. > //PetscPrintf(PETSC_COMM_WORLD,"%5i\t",closure[jj]); > indexXYZ = dim*(closure[jj]-vStart);//Linear index of x-coordinate in > the xyz_el array. > > *(x+vertices) = xyz_el[indexXYZ]; > *(y+vertices) = xyz_el[indexXYZ+1];//Extract Y-coordinates of the > current vertex. > *(z+vertices) = xyz_el[indexXYZ+2];//Extract Y-coordinates of the > current vertex. > *(EDOF + kk) = indexXYZ; > *(EDOF + kk+1) = indexXYZ+1; > *(EDOF + kk+2) = indexXYZ+2; > kk+=3; > vertices++;//Update vertex counter. > }//end if > else if(eStart <= closure[jj] && closure[jj]< eEnd){//Check for edge > ID's > edges++; > }//end else ifindexK > else if(fStart <= closure[jj] && closure[jj]< fEnd){//Check for face > ID's > faces++; > }//end else if > else if(cStart <= closure[jj] && closure[jj]< cEnd){//Check for cell > ID's > cells++; > }//end else if > }//end "jj" loop. > ierr = tetra4(x,y,z,&preKEtetra4,&matC,&KE); CHKERRQ(ierr); //Generate the > element stiffness matrix for this cell. > ierr = MatDenseGetArray(KE,&KEdata); CHKERRQ(ierr); > ierr = MatSetValues(K,12,EDOF,12,EDOF,KEdata,ADD_VALUES); > CHKERRQ(ierr);//WARNING: HARDCODED FOR TETRAHEDRAL P1 ELEMENTS ONLY > !!!!!!!!!!!!!!!!!!!!!!! > ierr = MatDenseRestoreArray(KE,&KEdata); CHKERRQ(ierr); > ierr = DMPlexRestoreTransitiveClosure(dm, ii,useCone, &size_closure, > &closure); CHKERRQ(ierr); > }//end "ii" loop. > ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > //ierr = MatView(K,PETSC_VIEWER_DRAW_WORLD); CHKERRQ(ierr); > //END GLOBAL STIFFNESS MATRIX > BUILDER=========================================== > t2 = MPI_Wtime(); > PetscPrintf(PETSC_COMM_WORLD,"K build time: %10f\n",t2-t1); > > > > > > > > > t1 = MPI_Wtime(); > //BEGIN BOUNDARY CONDITION > ENFORCEMENT========================================== > IS TrianglesIS, physicalsurfaceID;//, VerticesIS; > PetscInt numsurfvals, > //numRows, > dof_offset,numTri; > const PetscInt *surfvals, > //*pinZID, > *TriangleID; > PetscScalar diag =1; > PetscReal area,force; > //NOTE: Petsc can read/assign labels. Eeach label may posses multiple > "values." > //These values act as tags within a tag. > //IMPORTANT NOTE: The below line needs a safety. If a mesh that does not > feature > //face sets is imported, the code in its current state will crash!!!. This is > currently > //hardcoded for the test mesh. > ierr = DMGetLabel(dm, "Face Sets", &physicalgroups); CHKERRQ(ierr);//Inspects > Physical surface groups defined by gmsh (if any). > ierr = DMLabelGetValueIS(physicalgroups, &physicalsurfaceID); > CHKERRQ(ierr);//Gets the physical surface ID's defined in gmsh (as specified > in the .geo file). > ierr = ISGetIndices(physicalsurfaceID,&surfvals); CHKERRQ(ierr);//Get a > pointer to the actual surface values. > ierr = DMLabelGetNumValues(physicalgroups, &numsurfvals); > CHKERRQ(ierr);//Gets the number of different values that the label assigns. > for(ii=0;ii<numsurfvals;ii++){//loop through the values under the label. > //PetscPrintf(PETSC_COMM_WORLD,"Values = %5i\n",surfvals[ii]); > //PROBLEM: The surface values are hardcoded in the gmsh file. We need to > adopt standard "codes" > //that we can give to users when they make their meshes so that this code > recognizes the Type > // of boundary conditions that are to be imposed. > if(surfvals[ii] == pinXcode){ > dof_offset = 0; > dirichletBC = PETSC_TRUE; > }//end if. > else if(surfvals[ii] == pinZcode){ > dof_offset = 2; > dirichletBC = PETSC_TRUE; > }//end else if. > else if(surfvals[ii] == forceZcode){ > dof_offset = 2; > neumannBC = PETSC_TRUE; > }//end else if. > > ierr = DMLabelGetStratumIS(physicalgroups, surfvals[ii], &TrianglesIS); > CHKERRQ(ierr);//Get the ID's (as an IS) of the surfaces belonging to value 11. > //PROBLEM: DMPlexGetConeRecursiveVertices returns an array with repeated > node ID's. For each repetition, the lines that enforce BC's unnecessarily > re-run. > ierr = ISGetSize(TrianglesIS,&numTri); CHKERRQ(ierr); > ierr = ISGetIndices(TrianglesIS,&TriangleID); CHKERRQ(ierr);//Get a pointer > to the actual surface values. > for(kk=0;kk<numTri;kk++){ > ierr = DMPlexGetTransitiveClosure(dm, TriangleID[kk], useCone, > &size_closure, &closure); CHKERRQ(ierr); > if(neumannBC){ > ierr = DMPlexComputeCellGeometryFVM(dm, TriangleID[kk], > &area,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr); > force = traction*area/3;//WARNING: The 3 here is hardcoded for a purely > tetrahedral mesh only!!!!!!!!!! > } > for(jj=0;jj<(2*size_closure);jj+=2){ > //PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, > PetscReal *vol, PetscReal centroid[], PetscReal normal[]) > if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices. > indexK = dof*(closure[jj] - vStart) + dof_offset; //Compute the dof > ID's in the K matrix. > if(dirichletBC){//Boundary conditions requiring an edit of K matrix. > ierr = MatZeroRows(K,1,&indexK,diag,NULL,NULL); CHKERRQ(ierr); > }//end if. > else if(neumannBC){//Boundary conditions requiring an edit of RHS > vector. > ierr = VecSetValue(F,indexK,force,ADD_VALUES); CHKERRQ(ierr); > }// end else if. > }//end if. > }//end "jj" loop. > ierr = DMPlexRestoreTransitiveClosure(dm, closure[jj],useCone, > &size_closure, &closure); CHKERRQ(ierr); > }//end "kk" loop. > ierr = ISRestoreIndices(TrianglesIS,&TriangleID); CHKERRQ(ierr); > > /* > ierr = DMPlexGetConeRecursiveVertices(dm, TrianglesIS, &VerticesIS); > CHKERRQ(ierr);//Get the ID's (as an IS) of the vertices that make up the > surfaces of value 11. > ierr = ISGetSize(VerticesIS,&numRows); CHKERRQ(ierr);//Get number of > flagged vertices (this includes repeated indices for faces that share nodes). > ierr = ISGetIndices(VerticesIS,&pinZID); CHKERRQ(ierr);//Get a pointer to > the actual surface values. > if(dirichletBC){//Boundary conditions requiring an edit of K matrix. > for(kk=0;kk<numRows;kk++){ > indexK = 3*(pinZID[kk] - vStart) + dof_offset; //Compute the dof ID's > in the K matrix. (NOTE: the 3* ishardcoded for 3 degrees of freedom, tie this > to a variable in the FUTURE.) > ierr = MatZeroRows(K,1,&indexK,diag,NULL,NULL); CHKERRQ(ierr); > }//end "kk" loop. > }//end if. > else if(neumannBC){//Boundary conditions requiring an edit of RHS vector. > for(kk=0;kk<numRows;kk++){ > indexK = 3*(pinZID[kk] - vStart) + dof_offset; > ierr = VecSetValue(F,indexK,traction,INSERT_VALUES); CHKERRQ(ierr); > }//end "kk" loop. > }// end else if. > ierr = ISRestoreIndices(VerticesIS,&pinZID); CHKERRQ(ierr); > */ > dirichletBC = PETSC_FALSE; > neumannBC = PETSC_FALSE; > }//end "ii" loop. > ierr = ISRestoreIndices(physicalsurfaceID,&surfvals); CHKERRQ(ierr); > //ierr = ISRestoreIndices(VerticesIS,&pinZID); CHKERRQ(ierr); > ierr = ISDestroy(&physicalsurfaceID); CHKERRQ(ierr); > //ierr = ISDestroy(&VerticesIS); CHKERRQ(ierr); > ierr = ISDestroy(&TrianglesIS); CHKERRQ(ierr); > //END BOUNDARY CONDITION > ENFORCEMENT============================================ > t2 = MPI_Wtime(); > PetscPrintf(PETSC_COMM_WORLD,"BC imposition time: %10f\n",t2-t1); > > /* > PetscInt kk = 0; > for(ii=vStart;ii<vEnd;ii++){ > kk++; > PetscPrintf(PETSC_COMM_WORLD,"Vertex #%4i\t x = %10.9f\ty = %10.9f\tz = > %10.9f\n",ii,xyz_el[3*kk],xyz_el[3*kk+1],xyz_el[3*kk+2]); > }// end "ii" loop. > */ > > t1 = MPI_Wtime(); > //SOLVER======================================================================== > ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr); > ierr = KSPSetOperators(ksp,K,K); CHKERRQ(ierr); > ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr); > ierr = KSPSolve(ksp,F,U); CHKERRQ(ierr); > t2 = MPI_Wtime(); > //ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); > //SOLVER======================================================================== > t2 = MPI_Wtime(); > PetscPrintf(PETSC_COMM_WORLD,"Solver time: %10f\n",t2-t1); > ierr = VecRestoreArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's > data. > > //BEGIN MAX/MIN > DISPLACEMENTS=================================================== > IS ISux,ISuy,ISuz; > Vec UX,UY,UZ; > PetscReal UXmax,UYmax,UZmax,UXmin,UYmin,UZmin; > ierr = ISCreateStride(PETSC_COMM_WORLD,nv,0,3,&ISux); CHKERRQ(ierr); > ierr = ISCreateStride(PETSC_COMM_WORLD,nv,1,3,&ISuy); CHKERRQ(ierr); > ierr = ISCreateStride(PETSC_COMM_WORLD,nv,2,3,&ISuz); CHKERRQ(ierr); > > //PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y) > ierr = VecGetSubVector(U,ISux,&UX); CHKERRQ(ierr); > ierr = VecGetSubVector(U,ISuy,&UY); CHKERRQ(ierr); > ierr = VecGetSubVector(U,ISuz,&UZ); CHKERRQ(ierr); > > //PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val) > ierr = VecMax(UX,PETSC_NULL,&UXmax); CHKERRQ(ierr); > ierr = VecMax(UY,PETSC_NULL,&UYmax); CHKERRQ(ierr); > ierr = VecMax(UZ,PETSC_NULL,&UZmax); CHKERRQ(ierr); > > ierr = VecMin(UX,PETSC_NULL,&UXmin); CHKERRQ(ierr); > ierr = VecMin(UY,PETSC_NULL,&UYmin); CHKERRQ(ierr); > ierr = VecMin(UZ,PETSC_NULL,&UZmin); CHKERRQ(ierr); > > PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= ux <= %10f\n",UXmin,UXmax); > PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= uy <= %10f\n",UYmin,UYmax); > PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= uz <= %10f\n",UZmin,UZmax); > > > > > //BEGIN OUTPUT > SOLUTION========================================================= > if(saveASCII){ > ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"XYZ.txt",&XYZviewer); > ierr = VecView(XYZ,XYZviewer); CHKERRQ(ierr); > ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"U.txt",&XYZpUviewer); > ierr = VecView(U,XYZpUviewer); CHKERRQ(ierr); > PetscViewerDestroy(&XYZviewer); PetscViewerDestroy(&XYZpUviewer); > > }//end if. > if(saveVTK){ > const char *meshfile = "starting_mesh.vtk", > *deformedfile = "deformed_mesh.vtk"; > ierr = > PetscViewerVTKOpen(PETSC_COMM_WORLD,meshfile,FILE_MODE_WRITE,&XYZviewer); > CHKERRQ(ierr); > //PetscErrorCode DMSetAuxiliaryVec(DM dm, DMLabel label, PetscInt value, > Vec aux) > DMLabel UXlabel,UYlabel, UZlabel; > //PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel > *label) > ierr = DMLabelCreate(PETSC_COMM_WORLD, "X-Displacement", &UXlabel); > CHKERRQ(ierr); > ierr = DMLabelCreate(PETSC_COMM_WORLD, "Y-Displacement", &UYlabel); > CHKERRQ(ierr); > ierr = DMLabelCreate(PETSC_COMM_WORLD, "Z-Displacement", &UZlabel); > CHKERRQ(ierr); > ierr = DMSetAuxiliaryVec(dm,UXlabel, 1, UX); CHKERRQ(ierr); > ierr = DMSetAuxiliaryVec(dm,UYlabel, 1, UY); CHKERRQ(ierr); > ierr = DMSetAuxiliaryVec(dm,UZlabel, 1, UZ); CHKERRQ(ierr); > //PetscErrorCode PetscViewerVTKAddField(PetscViewer viewer,PetscObject > dm,PetscErrorCode > (*PetscViewerVTKWriteFunction)(PetscObject,PetscViewer),PetscInt > fieldnum,PetscViewerVTKFieldType fieldtype,PetscBool checkdm,PetscObject vec) > > > > //ierr = PetscViewerVTKAddField(XYZviewer, dm,PetscErrorCode > (*PetscViewerVTKWriteFunction)(Vec,PetscViewer),PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,UX); > ierr = PetscViewerVTKAddField(XYZviewer, > (PetscObject)dm,&PetscViewerVTKWriteFunction,PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,(PetscObject)UX); > > > ierr = DMPlexVTKWriteAll((PetscObject)dm, XYZviewer); CHKERRQ(ierr); > ierr = VecAXPY(XYZ,1,U); CHKERRQ(ierr);//Add displacement field to the mesh > coordinates to deform. > ierr = > PetscViewerVTKOpen(PETSC_COMM_WORLD,deformedfile,FILE_MODE_WRITE,&XYZpUviewer); > CHKERRQ(ierr); > ierr = DMPlexVTKWriteAll((PetscObject)dm, XYZpUviewer); CHKERRQ(ierr);// > PetscViewerDestroy(&XYZviewer); PetscViewerDestroy(&XYZpUviewer); > > }//end else if. > else{ > ierr = PetscPrintf(PETSC_COMM_WORLD,"No output format specified! Files not > saved.\n"); CHKERRQ(ierr); > }//end else. > > > //END OUTPUT > SOLUTION=========================================================== > VecDestroy(&UX); ISDestroy(&ISux); > VecDestroy(&UY); ISDestroy(&ISuy); > VecDestroy(&UZ); ISDestroy(&ISuz); > //END MAX/MIN > DISPLACEMENTS===================================================== > > > //CLEANUP===================================================================== > DMDestroy(&dm); > KSPDestroy(&ksp); > MatDestroy(&K); MatDestroy(&KE); MatDestroy(&matC); > //MatDestroy(preKEtetra4.matB); MatDestroy(preKEtetra4.matBTCB); > VecDestroy(&U); VecDestroy(&F); > > //DMLabelDestroy(&physicalgroups);//Destroyig the DM destroys the label. > > //CLEANUP===================================================================== > //PetscErrorCode PetscMallocDump(FILE *fp) > //ierr = PetscMallocDump(NULL); > return PetscFinalize();//And the machine shall rest.... > }//end main. > > PetscErrorCode tetra4(PetscScalar* X,PetscScalar* Y, PetscScalar* Z,struct > preKE *P, Mat* matC, Mat* KE){ > //INPUTS: > //X: Global X coordinates of the elemental nodes. > //Y: Global Y coordinates of the elemental nodes. > //Z: Global Z coordinates of the elemental nodes. > //J: Jacobian matrix. > //invJ: Inverse Jacobian matrix. > PetscErrorCode ierr; > //For current quadrature point, get dPsi/dXi_i Xi_i = {Xi,Eta,Zeta} > /* > P->dPsidXi[0] = +1.; P->dPsidEta[0] = 0.0; P->dPsidZeta[0] = 0.0; > P->dPsidXi[1] = 0.0; P->dPsidEta[1] = +1.; P->dPsidZeta[1] = 0.0; > P->dPsidXi[2] = 0.0; P->dPsidEta[2] = 0.0; P->dPsidZeta[2] = +1.; > P->dPsidXi[3] = -1.; P->dPsidEta[3] = -1.; P->dPsidZeta[3] = -1.; > */ > //Populate the Jacobian matrix. > P->J[0][0] = X[0] - X[3]; > P->J[0][1] = Y[0] - Y[3]; > P->J[0][2] = Z[0] - Z[3]; > P->J[1][0] = X[1] - X[3]; > P->J[1][1] = Y[1] - Y[3]; > P->J[1][2] = Z[1] - Z[3]; > P->J[2][0] = X[2] - X[3]; > P->J[2][1] = Y[2] - Y[3]; > P->J[2][2] = Z[2] - Z[3]; > > //Determinant of the 3x3 Jacobian. (Expansion along 1st row). > P->minor00 = P->J[1][1]*P->J[2][2] - P->J[2][1]*P->J[1][2];//Reuse when > finding InvJ. > P->minor01 = P->J[1][0]*P->J[2][2] - P->J[2][0]*P->J[1][2];//Reuse when > finding InvJ. > P->minor02 = P->J[1][0]*P->J[2][1] - P->J[2][0]*P->J[1][1];//Reuse when > finding InvJ. > P->detJ = P->J[0][0]*P->minor00 - P->J[0][1]*P->minor01 + > P->J[0][2]*P->minor02; > //Inverse of the 3x3 Jacobian > P->invJ[0][0] = +P->minor00/P->detJ;//Reuse precomputed minor. > P->invJ[0][1] = -(P->J[0][1]*P->J[2][2] - P->J[0][2]*P->J[2][1])/P->detJ; > P->invJ[0][2] = +(P->J[0][1]*P->J[1][2] - P->J[1][1]*P->J[0][2])/P->detJ; > P->invJ[1][0] = -P->minor01/P->detJ;//Reuse precomputed minor. > P->invJ[1][1] = +(P->J[0][0]*P->J[2][2] - P->J[0][2]*P->J[2][0])/P->detJ; > P->invJ[1][2] = -(P->J[0][0]*P->J[1][2] - P->J[1][0]*P->J[0][2])/P->detJ; > P->invJ[2][0] = +P->minor02/P->detJ;//Reuse precomputed minor. > P->invJ[2][1] = -(P->J[0][0]*P->J[2][1] - P->J[0][1]*P->J[2][0])/P->detJ; > P->invJ[2][2] = +(P->J[0][0]*P->J[1][1] - P->J[0][1]*P->J[1][0])/P->detJ; > > //*****************STRAIN MATRIX (B)************************************** > for(P->m=0;P->m<P->N;P->m++){//Scan all shape functions. > > P->x_in = 0 + P->m*3;//Every 3rd column starting at 0 > P->y_in = P->x_in +1;//Every 3rd column starting at 1 > P->z_in = P->y_in +1;//Every 3rd column starting at 2 > > P->dX[0] = P->invJ[0][0]*P->dPsidXi[P->m] + > P->invJ[0][1]*P->dPsidEta[P->m] + P->invJ[0][2]*P->dPsidZeta[P->m]; > P->dY[0] = P->invJ[1][0]*P->dPsidXi[P->m] + > P->invJ[1][1]*P->dPsidEta[P->m] + P->invJ[1][2]*P->dPsidZeta[P->m]; > P->dZ[0] = P->invJ[2][0]*P->dPsidXi[P->m] + > P->invJ[2][1]*P->dPsidEta[P->m] + P->invJ[2][2]*P->dPsidZeta[P->m]; > > P->dX[1] = P->dZ[0]; P->dX[2] = P->dY[0]; > P->dY[1] = P->dZ[0]; P->dY[2] = P->dX[0]; > P->dZ[1] = P->dX[0]; P->dZ[2] = P->dY[0]; > > ierr = > MatSetValues(P->matB,3,P->x_insert,1,&(P->x_in),P->dX,INSERT_VALUES); > CHKERRQ(ierr); > ierr = > MatSetValues(P->matB,3,P->y_insert,1,&(P->y_in),P->dY,INSERT_VALUES); > CHKERRQ(ierr); > ierr = > MatSetValues(P->matB,3,P->z_insert,1,&(P->z_in),P->dZ,INSERT_VALUES); > CHKERRQ(ierr); > > }//end "m" loop. > ierr = MatAssemblyBegin(P->matB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > ierr = MatAssemblyEnd(P->matB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > //*****************STRAIN MATRIX (B)************************************** > > //Compute the matrix product B^t*C*B, scale it by the quadrature > weights and add to KE. > P->weight = -P->detJ/6; > > ierr = MatZeroEntries(*KE); CHKERRQ(ierr); > ierr = > MatPtAP(*matC,P->matB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(P->matBTCB));CHKERRQ(ierr); > ierr = MatScale(P->matBTCB,P->weight); CHKERRQ(ierr); > ierr = MatAssemblyBegin(P->matBTCB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > ierr = MatAssemblyEnd(P->matBTCB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > ierr = MatAXPY(*KE,1,P->matBTCB,DIFFERENT_NONZERO_PATTERN); > CHKERRQ(ierr);//Add contribution of current quadrature point to KE. > > //ierr = > MatPtAP(*matC,P->matB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,KE);CHKERRQ(ierr); > //ierr = MatScale(*KE,P->weight); CHKERRQ(ierr); > > ierr = MatAssemblyBegin(*KE,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > ierr = MatAssemblyEnd(*KE,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > > //Cleanup > return ierr; > }//end tetra4. > > PetscErrorCode ConstitutiveMatrix(Mat *matC,const char* type,PetscInt > materialID){ > PetscErrorCode ierr; > PetscBool isotropic = PETSC_FALSE, > orthotropic = PETSC_FALSE; > //PetscErrorCode PetscStrcmp(const char a[],const char b[],PetscBool *flg) > ierr = PetscStrcmp(type,"isotropic",&isotropic); > ierr = PetscStrcmp(type,"orthotropic",&orthotropic); > ierr = MatCreate(PETSC_COMM_WORLD,matC); CHKERRQ(ierr); > ierr = MatSetSizes(*matC,PETSC_DECIDE,PETSC_DECIDE,6,6); CHKERRQ(ierr); > ierr = MatSetType(*matC,MATAIJ); CHKERRQ(ierr); > ierr = MatSetUp(*matC); CHKERRQ(ierr); > > if(isotropic){ > PetscReal E,nu, M,L,vals[3]; > switch(materialID){ > case 0://Hardcoded properties for isotropic material #0 > E = 200; > nu = 1./3; > break; > case 1://Hardcoded properties for isotropic material #1 > E = 96; > nu = 1./3; > break; > }//end switch. > M = E/(2*(1+nu)),//Lame's constant 1 ("mu"). > L = E*nu/((1+nu)*(1-2*nu));//Lame's constant 2 ("lambda"). > //PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt > idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv) > PetscInt idxn[3] = {0,1,2}; > vals[0] = L+2*M; vals[1] = L; vals[2] = vals[1]; > ierr = MatSetValues(*matC,1,&idxn[0],3,idxn,vals,INSERT_VALUES); > CHKERRQ(ierr); > vals[1] = vals[0]; vals[0] = vals[2]; > ierr = MatSetValues(*matC,1,&idxn[1],3,idxn,vals,INSERT_VALUES); > CHKERRQ(ierr); > vals[2] = vals[1]; vals[1] = vals[0]; > ierr = MatSetValues(*matC,1,&idxn[2],3,idxn,vals,INSERT_VALUES); > CHKERRQ(ierr); > ierr = MatSetValue(*matC,3,3,M,INSERT_VALUES); CHKERRQ(ierr); > ierr = MatSetValue(*matC,4,4,M,INSERT_VALUES); CHKERRQ(ierr); > ierr = MatSetValue(*matC,5,5,M,INSERT_VALUES); CHKERRQ(ierr); > }//end if. > /* > else if(orthotropic){ > switch(materialID){ > case 0: > break; > case 1: > break; > }//end switch. > }//end else if. > */ > ierr = MatAssemblyBegin(*matC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > ierr = MatAssemblyEnd(*matC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); > //MatView(*matC,0); > return ierr; > }//End ConstitutiveMatrix > > PetscErrorCode InitializeKEpreallocation(struct preKE *P,const char* type){ > PetscErrorCode ierr; > PetscBool istetra4 = PETSC_FALSE, > ishex8 = PETSC_FALSE; > ierr = PetscStrcmp(type,"tetra4",&istetra4); CHKERRQ(ierr); > ierr = PetscStrcmp(type,"hex8",&ishex8); CHKERRQ(ierr); > if(istetra4){ > P->sizeKE = 12; > P->N = 4; > }//end if. > else if(ishex8){ > P->sizeKE = 24; > P->N = 8; > }//end else if. > > > P->x_insert[0] = 0; P->x_insert[1] = 3; P->x_insert[2] = 5; > P->y_insert[0] = 1; P->y_insert[1] = 4; P->y_insert[2] = 5; > P->z_insert[0] = 2; P->z_insert[1] = 3; P->z_insert[2] = 4; > //Allocate memory for the differentiated shape function vectors. > ierr = PetscMalloc1(P->N,&(P->dPsidXi)); CHKERRQ(ierr); > ierr = PetscMalloc1(P->N,&(P->dPsidEta)); CHKERRQ(ierr); > ierr = PetscMalloc1(P->N,&(P->dPsidZeta)); CHKERRQ(ierr); > > P->dPsidXi[0] = +1.; P->dPsidEta[0] = 0.0; P->dPsidZeta[0] = 0.0; > P->dPsidXi[1] = 0.0; P->dPsidEta[1] = +1.; P->dPsidZeta[1] = 0.0; > P->dPsidXi[2] = 0.0; P->dPsidEta[2] = 0.0; P->dPsidZeta[2] = +1.; > P->dPsidXi[3] = -1.; P->dPsidEta[3] = -1.; P->dPsidZeta[3] = -1.; > > > //Strain matrix. > ierr = MatCreate(PETSC_COMM_WORLD,&(P->matB)); CHKERRQ(ierr); > ierr = MatSetSizes(P->matB,PETSC_DECIDE,PETSC_DECIDE,6,P->sizeKE); > CHKERRQ(ierr);//Hardcoded > ierr = MatSetType(P->matB,MATAIJ); CHKERRQ(ierr); > ierr = MatSetUp(P->matB); CHKERRQ(ierr); > > //Contribution matrix. > ierr = MatCreate(PETSC_COMM_WORLD,&(P->matBTCB)); CHKERRQ(ierr); > ierr = > MatSetSizes(P->matBTCB,PETSC_DECIDE,PETSC_DECIDE,P->sizeKE,P->sizeKE); > CHKERRQ(ierr); > ierr = MatSetType(P->matBTCB,MATAIJ); CHKERRQ(ierr); > ierr = MatSetUp(P->matBTCB); CHKERRQ(ierr); > > //Element stiffness matrix. > //ierr = MatCreateSeqDense(PETSC_COMM_SELF,12,12,NULL,&KE); CHKERRQ(ierr); > //PARALLEL > > return ierr; > }