On Thu, Jan 6, 2022 at 5:15 PM Ferrand, Jesus A. <ferra...@my.erau.edu> wrote:
> Jed: > > DMPlexLabelComplete() has allowed me to speed up my code significantly > (Many thanks!). > > I did not use DMAddBoundary() though. > I figured I could obtain Index Sets (IS's) from the DAG for different > depths and then IS's for the points that were flagged in Gmsh (after > calling DMPlexLabelComplete()). > I then intersected both IS's using ISIntersect() to get the DAG points > corresponding to just vertices (and flagged by Gmsh) for Dirichlet BC's, > and DAG points that are Faces and flagged by Gmsh for Neumann BC's. I then > use the intersected IS to edit a Mat and a RHS Vec manually. I did further > profiling and have found the PetsSections are now the next biggest > overhead. > > For Dirichlet BC's I make an array of vertex ID's and call > MatSetZeroRows() to impose BC's on them through my K matrix. And yes, I'm > solving the elasticity PDE. For Neumann BC's I use > DMPlexGetRecursiveVertices() to edit my RHS vector. > I cannot find a function named DMPlexGetRecursiveVertices(). > I want to keep the PetscSections since they preallocate my matrix rather > well (the one from DMCreateMatrix()) but at the same time would like to > remove them since they add overhead. Do you think DMAddboundary() with the > function call will be faster than my single calls to MatSetZeroRows() and > DMPlexGetRecursiveVertices() ? > PetscSection is really simple. Are you sure you are measuring long times there? What are you using it to do? Thanks, Matt > ------------------------------ > *From:* Jed Brown <j...@jedbrown.org> > *Sent:* Wednesday, January 5, 2022 5:44 PM > *To:* Ferrand, Jesus A. <ferra...@my.erau.edu> > *Cc:* petsc-users <petsc-users@mcs.anl.gov> > *Subject:* Re: [EXTERNAL] Re: [petsc-users] DM misuse causes massive > memory leak? > > For something like displacement (and this sounds like elasticity), I would > recommend using one field with three components. You can constrain a subset > of the components to implement slip conditions. > > You can use DMPlexLabelComplete(dm, label) to propagate those face labels > to vertices. > > "Ferrand, Jesus A." <ferra...@my.erau.edu> writes: > > > Thanks for the reply (I hit reply all this time). > > > > So, I set 3 fields: > > /* > > ierr = PetscSectionSetNumFields(s,dof); CHKERRQ(ierr); > > ierr = PetscSectionSetFieldName(s,0, "X-Displacement"); CHKERRQ(ierr); > //Field ID is 0 > > ierr = PetscSectionSetFieldName(s,1, "Y-Displacement"); CHKERRQ(ierr); > //Field ID is 1 > > ierr = PetscSectionSetFieldName(s,2, "Z-Displacement"); CHKERRQ(ierr); > //Field ID is 2 > > */ > > > > I then loop through the vertices of my DMPlex > > > > /* > > for(ii = vStart; ii < vEnd; ii++){//Vertex loop. > > ierr = PetscSectionSetDof(s, ii, dof); CHKERRQ(ierr); > > ierr = PetscSectionSetFieldDof(s,ii,0,1); CHKERRQ(ierr);//One > X-displacement per vertex (1 dof) > > ierr = PetscSectionSetFieldDof(s,ii,1,1); CHKERRQ(ierr);//One > Y-displacement per vertex (1 dof) > > ierr = PetscSectionSetFieldDof(s,ii,2,1); CHKERRQ(ierr);//One > Z-displacement per vertex (1 dof) > > }//Sets x, y, and z displacements as dofs. > > */ > > > > I only associated fields with vertices, not with any other points in the > DAG. Regarding the use of DMAddBoundary(), I mostly copied the usage shown > in SNES example 77. I modified the function definition to simply set the > dof to 0.0 as opposed to the coordinates. Below "physicalgroups" is the > DMLabel that I got from gmsh, this flags Face points, not vertices. That is > why I think the error log suggests that fields were never set. > > > > ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", physicalgroups, 1, > &surfvals[ii], fieldID, 0, NULL, (void (*)(void)) coordinates, NULL, NULL, > NULL); CHKERRQ(ierr); > > PetscErrorCode coordinates(PetscInt dim, PetscReal time, const PetscReal > x[], PetscInt Nf, PetscScalar *u, void *ctx){ > > const PetscInt Ncomp = dim; > > PetscInt comp; > > for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0; > > return 0; > > } > > > > > > ________________________________ > > From: Jed Brown <j...@jedbrown.org> > > Sent: Wednesday, January 5, 2022 12:36 AM > > To: Ferrand, Jesus A. <ferra...@my.erau.edu> > > Cc: petsc-users <petsc-users@mcs.anl.gov> > > Subject: Re: [EXTERNAL] Re: [petsc-users] DM misuse causes massive > memory leak? > > > > Please "reply all" include the list in the future. > > > > "Ferrand, Jesus A." <ferra...@my.erau.edu> writes: > > > >> Forgot to say thanks for the reply (my bad). > >> Yes, I was indeed forgetting to pre-allocate the sparse matrices when > doing them by hand (complacency seeded by MATLAB's zeros()). Thank you, > Jed, and Jeremy, for the hints! > >> > >> I have more questions, these ones about boundary conditions (I think > these are for Matt). > >> In my current code I set Dirichlet conditions directly on a Mat by > calling MatSetZeroRows(). I profiled my code and found the part that > applies them to be unnacceptably slow. In response, I've been trying to > better pre-allocate Mats using PetscSections. I have found documentation > for PetscSectionSetDof(), PetscSectionSetNumFields(), > PetscSectionSetFieldName(), and PetscSectionSetFieldDof(), to set the size > of my Mats and Vecs by calling DMSetLocalSection() followed by > DMCreateMatrix() and DMCreateLocalVector() to get a RHS vector. This seems > faster. > >> > >> In PetscSection, what is the difference between a "field" and a > "component"? For example, I could have one field "Velocity" with three > components ux, uy, and uz or perhaps three fields ux, uy, and uz each with > a default component? > > > > It's just how you name them and how they appear in output. Usually > "velocity" is better as a field with three components, but fields with > other meaning (and perhaps different finite element spaces), such as > pressure, would be different fields. Different components are always in the > same FE space. > > > >> I am struggling now to impose boundary conditions after constraining > dofs using PetscSection. My understanding is that constraining dof's > reduces the size of the DM's matrix but it does not give the DM knowledge > of what values the constrained dofs should have, right? > >> > >> I know that there is DMAddBoundary(), but I am unsure of how to use it. > From Gmsh I have a mesh with surface boundaries flagged. I'm not sure > whether DMAddBoundary()will constrain the face, edge, or vertex points when > I give it the DMLabel from Gmsh. (I need specific dof on the vertices to be > constrained). I did some testing and I think DMAddBoundary() attempts to > constrain the Face points (see error log below). I only associated fields > with the vertices but not the Faces. I can extract the vertex points from > the face label using DMPlexGetConeRecursiveVertices() but the output IS has > repeated entries for the vertex points (many faces share the same vertex). > Is there an easier way to get the vertex points from a gmsh surface tag? > > > > How did you call DMAddBoundary()? Are you using DM_BC_ESSENTIAL and a > callback function that provides the inhomogeneous boundary condition? > > > >> I'm sorry this is a mouthful. > >> > >> [0]PETSC ERROR: --------------------- Error Message > -------------------------------------------------------------- > >> [0]PETSC ERROR: Argument out of range > >> [0]PETSC ERROR: Field number 2 must be in [0, 0) > > > > It looks like you haven't added these fields yet. > > > >> [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble > shooting. > >> [0]PETSC ERROR: Petsc Release Version 3.16.0, unknown > >> [0]PETSC ERROR: ./gmsh4 on a linux-c-dbg named F86 by jesus Tue Jan 4 > 15:19:57 2022 > >> [0]PETSC ERROR: Configure options --with-32bits-pci-domain=1 > --with-debugging =1 > >> [0]PETSC ERROR: #1 DMGetField() at > /home/jesus/SAND/PETSc_install/petsc/src/dm/interface/dm.c:4803 > >> [0]PETSC ERROR: #2 DMCompleteBoundaryLabel_Internal() at > /home/jesus/SAND/PETSc_install/petsc/src/dm/interface/dm.c:5140 > >> [0]PETSC ERROR: #3 DMAddBoundary() at > /home/jesus/SAND/PETSc_install/petsc/src/dm/interface/dm.c:8561 > >> [0]PETSC ERROR: #4 main() at /home/jesus/SAND/FEA/3D/gmshBACKUP4.c:215 > >> [0]PETSC ERROR: No PETSc Option Table entries > >> [0]PETSC ERROR: ----------------End of Error Message -------send entire > error message to petsc-ma...@mcs.anl.gov---------- > >> > >> > >> > >> > >> > >> > >> ________________________________ > >> From: Jed Brown <j...@jedbrown.org> > >> Sent: Wednesday, December 29, 2021 5:55 PM > >> To: Ferrand, Jesus A. <ferra...@my.erau.edu>; petsc-users@mcs.anl.gov < > petsc-users@mcs.anl.gov> > >> Subject: [EXTERNAL] Re: [petsc-users] DM misuse causes massive memory > leak? > >> > >> CAUTION: This email originated outside of Embry-Riddle Aeronautical > University. Do not click links or open attachments unless you recognize the > sender and know the content is safe. > >> > >> > >> "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; > >>> } > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>