[petsc-users] [EXTERNAL] Re: DM misuse causes massive memory leak?
Matthew Knepley
knepley at gmail.com
Fri Jan 7 18:23:02 CST 2022
On Fri, Jan 7, 2022 at 7:10 PM Ferrand, Jesus A. <FERRANJ2 at my.erau.edu>
wrote:
> Matt:
>
> Thanks again for the prompt reply : )
>
> You were right regarding the PetscSections. I reprofiled my code where I
> use them and realized that I had DMCreateMatrix() in the mix. The brunt of
> the overhead came from it, not the PetscSections.
>
> I looked for ways to improve the time on DMCreateMatrix() and figured I
> could call DMSetMatType() to work with a BAIJ + MatSetValuesBlocked()
> instead of AIJ +MatSetValues() to leverage my dof blocks (ux, uy, uz)
> during assembly. My matrix assembly times have improved but DMCreateMatrix
> seems to take just as long as before. I noticed after calling
> DMGetBlockSize() that my DM is working with blocksize of 1. Weirdly enough,
> my call to KSPView() returns a block size of 3 (as it should be) for the
> matrix that came from DMCreateMatrix(). I further noticed after looking at
> the DM function index that there is no "DMSetBlockSize()," however there is
> in fact a DMGetBlockSize(). Anyways, any suggestions on how to speed up
> DMCreateMatrix()? I have tried calling both DMSetMatrixPreallocateOnly()
> and DMSetAdjacency() (with the flags for FEM) but barely got a speedup.
>
THe time here comes from traversing the mesh to see what needs to be
allocated. It is insensitive to block size. I know of no way to speed up
this process
without strong additional assumptions.
> Regarding DMPlexGetConeRecursiveVertices() vs.
> DMPlexGetTransitiveClosure() for BC's. I will (for now) side with the
> former as I am getting better times with it for the BC's at least. My gripe
> with DMPlexGetTransitiveClosure() is that it returns the whole closure as
> oppossed to subsets of it (i.e., for a cell, it returns the constituetn
> faces, edges, and vertices as oppossed to say, just vertices). Currently,
> I'm working with P1 elements only so I need only vertices. To me,
> therefore, the additional info seems wasteful for P1. For Pn (n>1)
> DMPlexGetTransitiveClosure() seems like the superior option though.
>
GetConeRecursive _is_ getting the entire closure and throwing away the
parts it does not want. It is exactly the same code with worse memory
management.
Thanks,
Matt
> ------------------------------
> *From:* Matthew Knepley <knepley at gmail.com>
> *Sent:* Thursday, January 6, 2022 5:39 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: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/>
>
--
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/20220107/bf3b7b37/attachment-0001.html>
More information about the petsc-users
mailing list