[petsc-users] [EXTERNAL] Re: DM misuse causes massive memory leak?
Matthew Knepley
knepley at gmail.com
Thu Jan 6 16:39:10 CST 2022
On Thu, Jan 6, 2022 at 5:30 PM Ferrand, Jesus A. <FERRANJ2 at my.erau.edu>
wrote:
> Matt:
>
> Thanks for the immediate reply.
> My apologies, I misspelled the function name, it should have been "
> DMPlexGetConeRecursiveVertices()."
>
Oh, this is a debugging function that Vaclav wrote. It should never be used
in production. It allocates memory all over the place.
If you want vertices, just get the closure and filter out the vertices:
PetscInt *closure = NULL;
PetscInt Ncl, Nclv;
DMPlexGetDepthStrautm(dm, &vStart, &vEnd);
DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &Ncl, &closure);
Nclv = 0;
for (cl = 0; cl < 2*Ncl; cl += 2) {
const PetscInt point = closure[cl];
if ((point >= vStart) && (point < vEnd)) closure[Nclv++] = point;
}
<closure[] contains Nclv vertices>
DMPlexRestoreTransitiveClosure(...);
> Regarding the PetscSection use: I am looping through every single point in
> the DAG of my mesh. For each point I am assigning dof using
> PetscSectionSetDof(). I am also assigning dof to the corresponding fields
> using PetscSectionSetFieldDof(). I took Jed's advice and made a single
> field with 3 components, I named all of them. So, I used
> PetscSectionSetNumFields(), PetscSectionSetFieldComponents(),
> PetscSectionSetFieldName(), and PetscSectionSetComponentName(). Finally, I
> proceed to PetscSectionSetUp(), and DMSetLocalSection(). In my timed code
> blocks I am including DMPlexGetDepthStratum(), and DMGetStratumSize(), and
> DMPlexGetChart() because I need their output to assign the dof to the
> PetscSection.
>
This should all take almost no time. There are no expensive operations
there.
Thanks,
Matt
> For what is worth, I have my configure set to --with-debugging = 1.
> ------------------------------
> *From:* Matthew Knepley <knepley at gmail.com>
> *Sent:* Thursday, January 6, 2022 5:20 PM
> *To:* Ferrand, Jesus A. <FERRANJ2 at my.erau.edu>
> *Cc:* Jed Brown <jed at jedbrown.org>; petsc-users <petsc-users at mcs.anl.gov>
> *Subject:* Re: [petsc-users] [EXTERNAL] Re: DM misuse causes massive
> memory leak?
>
> On Thu, Jan 6, 2022 at 5:15 PM Ferrand, Jesus A. <FERRANJ2 at 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 <jed at jedbrown.org>
> *Sent:* Wednesday, January 5, 2022 5:44 PM
> *To:* Ferrand, Jesus A. <FERRANJ2 at my.erau.edu>
> *Cc:* petsc-users <petsc-users at 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." <FERRANJ2 at 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 <jed at jedbrown.org>
> > Sent: Wednesday, January 5, 2022 12:36 AM
> > To: Ferrand, Jesus A. <FERRANJ2 at my.erau.edu>
> > Cc: petsc-users <petsc-users at 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." <FERRANJ2 at 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-maint at mcs.anl.gov----------
> >>
> >>
> >>
> >>
> >>
> >>
> >> ________________________________
> >> From: Jed Brown <jed at jedbrown.org>
> >> Sent: Wednesday, December 29, 2021 5:55 PM
> >> To: Ferrand, Jesus A. <FERRANJ2 at my.erau.edu>; petsc-users at mcs.anl.gov <
> petsc-users at 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." <FERRANJ2 at 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): ferranj2 at my.erau.edu
> >>>
> >>> jesus.ferrand at 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/>
>
--
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220106/19dd9966/attachment-0001.html>
More information about the petsc-users
mailing list