<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body dir="ltr">
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
Jed:<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
<br>
DMPlexLabelComplete() has allowed me to speed up my code significantly (Many thanks!).</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
I did not use DMAddBoundary() though.<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
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()).<br>
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.
<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
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.<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
<br>
</div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0); background-color: rgb(255, 255, 255);">
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() ?<br>
</div>
<div>
<div id="appendonsend"></div>
<div style="font-family: "Times New Roman", Times, serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
<hr tabindex="-1" style="display:inline-block; width:98%">
<div id="divRplyFwdMsg" dir="ltr"><font style="font-size: 11pt;" face="Calibri, sans-serif" color="#000000"><b>From:</b> Jed Brown <jed@jedbrown.org><br>
<b>Sent:</b> Wednesday, January 5, 2022 5:44 PM<br>
<b>To:</b> Ferrand, Jesus A. <FERRANJ2@my.erau.edu><br>
<b>Cc:</b> petsc-users <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [EXTERNAL] Re: [petsc-users] DM misuse causes massive memory leak?</font>
<div> </div>
</div>
<div class="BodyFragment"><font size="2"><span style="font-size:11pt">
<div class="PlainText">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.<br>
<br>
You can use DMPlexLabelComplete(dm, label) to propagate those face labels to vertices.<br>
<br>
"Ferrand, Jesus A." <FERRANJ2@my.erau.edu> writes:<br>
<br>
> Thanks for the reply (I hit reply all this time).<br>
><br>
> So, I set 3 fields:<br>
> /*<br>
> ierr = PetscSectionSetNumFields(s,dof); CHKERRQ(ierr);<br>
> ierr = PetscSectionSetFieldName(s,0, "X-Displacement"); CHKERRQ(ierr); //Field ID is 0<br>
> ierr = PetscSectionSetFieldName(s,1, "Y-Displacement"); CHKERRQ(ierr); //Field ID is 1<br>
> ierr = PetscSectionSetFieldName(s,2, "Z-Displacement"); CHKERRQ(ierr); //Field ID is 2<br>
> */<br>
><br>
> I then loop through the vertices of my DMPlex<br>
><br>
> /*<br>
> for(ii = vStart; ii < vEnd; ii++){//Vertex loop.<br>
> ierr = PetscSectionSetDof(s, ii, dof); CHKERRQ(ierr);<br>
> ierr = PetscSectionSetFieldDof(s,ii,0,1); CHKERRQ(ierr);//One X-displacement per vertex (1 dof)<br>
> ierr = PetscSectionSetFieldDof(s,ii,1,1); CHKERRQ(ierr);//One Y-displacement per vertex (1 dof)<br>
> ierr = PetscSectionSetFieldDof(s,ii,2,1); CHKERRQ(ierr);//One Z-displacement per vertex (1 dof)<br>
> }//Sets x, y, and z displacements as dofs.<br>
> */<br>
><br>
> 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.<br>
><br>
> ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "fixed", physicalgroups, 1, &surfvals[ii], fieldID, 0, NULL, (void (*)(void)) coordinates, NULL, NULL, NULL); CHKERRQ(ierr);<br>
> PetscErrorCode coordinates(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){<br>
> const PetscInt Ncomp = dim;<br>
> PetscInt comp;<br>
> for (comp = 0; comp < Ncomp; ++comp) u[comp] = 0.0;<br>
> return 0;<br>
> }<br>
><br>
><br>
> ________________________________<br>
> From: Jed Brown <jed@jedbrown.org><br>
> Sent: Wednesday, January 5, 2022 12:36 AM<br>
> To: Ferrand, Jesus A. <FERRANJ2@my.erau.edu><br>
> Cc: petsc-users <petsc-users@mcs.anl.gov><br>
> Subject: Re: [EXTERNAL] Re: [petsc-users] DM misuse causes massive memory leak?<br>
><br>
> Please "reply all" include the list in the future.<br>
><br>
> "Ferrand, Jesus A." <FERRANJ2@my.erau.edu> writes:<br>
><br>
>> Forgot to say thanks for the reply (my bad).<br>
>> 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!<br>
>><br>
>> I have more questions, these ones about boundary conditions (I think these are for Matt).<br>
>> 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.<br>
>><br>
>> 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?<br>
><br>
> 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.<br>
><br>
>> 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?<br>
>><br>
>> 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?<br>
><br>
> How did you call DMAddBoundary()? Are you using DM_BC_ESSENTIAL and a callback function that provides the inhomogeneous boundary condition?<br>
><br>
>> I'm sorry this is a mouthful.<br>
>><br>
>> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
>> [0]PETSC ERROR: Argument out of range<br>
>> [0]PETSC ERROR: Field number 2 must be in [0, 0)<br>
><br>
> It looks like you haven't added these fields yet.<br>
><br>
>> [0]PETSC ERROR: See <a href="https://petsc.org/release/faq/">https://petsc.org/release/faq/</a> for trouble shooting.<br>
>> [0]PETSC ERROR: Petsc Release Version 3.16.0, unknown<br>
>> [0]PETSC ERROR: ./gmsh4 on a linux-c-dbg named F86 by jesus Tue Jan 4 15:19:57 2022<br>
>> [0]PETSC ERROR: Configure options --with-32bits-pci-domain=1 --with-debugging =1<br>
>> [0]PETSC ERROR: #1 DMGetField() at /home/jesus/SAND/PETSc_install/petsc/src/dm/interface/dm.c:4803<br>
>> [0]PETSC ERROR: #2 DMCompleteBoundaryLabel_Internal() at /home/jesus/SAND/PETSc_install/petsc/src/dm/interface/dm.c:5140<br>
>> [0]PETSC ERROR: #3 DMAddBoundary() at /home/jesus/SAND/PETSc_install/petsc/src/dm/interface/dm.c:8561<br>
>> [0]PETSC ERROR: #4 main() at /home/jesus/SAND/FEA/3D/gmshBACKUP4.c:215<br>
>> [0]PETSC ERROR: No PETSc Option Table entries<br>
>> [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov----------<br>
>><br>
>><br>
>><br>
>><br>
>><br>
>><br>
>> ________________________________<br>
>> From: Jed Brown <jed@jedbrown.org><br>
>> Sent: Wednesday, December 29, 2021 5:55 PM<br>
>> To: Ferrand, Jesus A. <FERRANJ2@my.erau.edu>; petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov><br>
>> Subject: [EXTERNAL] Re: [petsc-users] DM misuse causes massive memory leak?<br>
>><br>
>> 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.<br>
>><br>
>><br>
>> "Ferrand, Jesus A." <FERRANJ2@my.erau.edu> writes:<br>
>><br>
>>> Dear PETSc Team:<br>
>>><br>
>>> 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().<br>
>><br>
>> 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.<br>
>><br>
>> <a href="https://petsc.org/release/docs/manual/performance/#memory-allocation-for-sparse-matrix-assembly">
https://petsc.org/release/docs/manual/performance/#memory-allocation-for-sparse-matrix-assembly</a><br>
>> <a href="https://petsc.org/release/docs/manual/mat/#sec-matsparse">https://petsc.org/release/docs/manual/mat/#sec-matsparse</a><br>
>><br>
>>> 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.<br>
>>><br>
>>> Machine Type: HP Laptop<br>
>>> C-compiler: Gnu C<br>
>>> OS: Ubuntu 20.04<br>
>>> PETSc version: 3.16.0<br>
>>> MPI Implementation: MPICH<br>
>>><br>
>>> Hope you all had a Merry Christmas and that you will have a happy and productive New Year. :D<br>
>>><br>
>>><br>
>>> Sincerely:<br>
>>><br>
>>> J.A. Ferrand<br>
>>><br>
>>> Embry-Riddle Aeronautical University - Daytona Beach FL<br>
>>><br>
>>> M.Sc. Aerospace Engineering | May 2022<br>
>>><br>
>>> B.Sc. Aerospace Engineering<br>
>>><br>
>>> B.Sc. Computational Mathematics<br>
>>><br>
>>><br>
>>><br>
>>> Sigma Gamma Tau<br>
>>><br>
>>> Tau Beta Pi<br>
>>><br>
>>> Honors Program<br>
>>><br>
>>><br>
>>><br>
>>> Phone: (386)-843-1829<br>
>>><br>
>>> Email(s): ferranj2@my.erau.edu<br>
>>><br>
>>> jesus.ferrand@gmail.com<br>
>>> //REFERENCE: <a href="https://github.com/FreeFem/FreeFem-sources/blob/master/plugin/mpi/PETSc-code.hpp">
https://github.com/FreeFem/FreeFem-sources/blob/master/plugin/mpi/PETSc-code.hpp</a><br>
>>> #include <petsc.h><br>
>>> static char help[] = "Imports a Gmsh mesh with boundary conditions and solves the elasticity equation.\n"<br>
>>> "Option prefix = opt_.\n";<br>
>>><br>
>>> struct preKE{//Preallocation before computing KE<br>
>>> Mat matB,<br>
>>> matBTCB;<br>
>>> //matKE;<br>
>>> PetscInt x_insert[3],<br>
>>> y_insert[3],<br>
>>> z_insert[3],<br>
>>> m,//Looping variables.<br>
>>> sizeKE,//size of the element stiffness matrix.<br>
>>> N,//Number of nodes in element.<br>
>>> x_in,y_in,z_in; //LI to index B matrix.<br>
>>> PetscReal J[3][3],//Jacobian matrix.<br>
>>> invJ[3][3],//Inverse of the Jacobian matrix.<br>
>>> detJ,//Determinant of the Jacobian.<br>
>>> dX[3],<br>
>>> dY[3],<br>
>>> dZ[3],<br>
>>> minor00,<br>
>>> minor01,<br>
>>> minor02,//Determinants of minors in a 3x3 matrix.<br>
>>> dPsidX, dPsidY, dPsidZ,//Shape function derivatives w.r.t global coordinates.<br>
>>> weight,//Multiplier of quadrature weights.<br>
>>> *dPsidXi,//Derivatives of shape functions w.r.t. Xi.<br>
>>> *dPsidEta,//Derivatives of shape functions w.r.t. Eta.<br>
>>> *dPsidZeta;//Derivatives of shape functions w.r.t Zeta.<br>
>>> PetscErrorCode ierr;<br>
>>> };//end struct.<br>
>>><br>
>>> //Function declarations.<br>
>>> extern PetscErrorCode tetra4(PetscScalar*, PetscScalar*, PetscScalar*,struct preKE*, Mat*, Mat*);<br>
>>> extern PetscErrorCode ConstitutiveMatrix(Mat*,const char*,PetscInt);<br>
>>> extern PetscErrorCode InitializeKEpreallocation(struct preKE*,const char*);<br>
>>><br>
>>> PetscErrorCode PetscViewerVTKWriteFunction(PetscObject vec,PetscViewer viewer){<br>
>>> PetscErrorCode ierr;<br>
>>> ierr = VecView((Vec)vec,viewer); CHKERRQ(ierr);<br>
>>> return ierr;<br>
>>> }<br>
>>><br>
>>><br>
>>><br>
>>><br>
>>> int main(int argc, char **args){<br>
>>> //DEFINITIONS OF PETSC's DMPLEX LINGO:<br>
>>> //POINT: A topology element (cell, face, edge, or vertex).<br>
>>> //CHART: It an interval from 0 to the number of "points." (the range of admissible linear indices)<br>
>>> //STRATUM: A subset of the "chart" which corresponds to all "points" at a given "level."<br>
>>> //LEVEL: This is either a "depth" or a "height".<br>
>>> //HEIGHT: Dimensionality of an element measured from 0D to 3D. Heights: cell = 0, face = 1, edge = 2, vertex = 3.<br>
>>> //DEPTH: Dimensionality of an element measured from 3D to 0D. Depths: cell = 3, face = 2, edge = 1, vertex = 0;<br>
>>> //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.<br>
>>> //STAR:<br>
>>> //STANDARD LABELS: These are default tags that DMPlex has for its topology. ("depth")<br>
>>> PetscErrorCode ierr;//Error tracking variable.<br>
>>> DM dm;//Distributed memory object (useful for managing grids.)<br>
>>> DMLabel physicalgroups;//Identifies user-specified tags in gmsh (to impose BC's).<br>
>>> DMPolytopeType celltype;//When looping through cells, determines its type (tetrahedron, pyramid, hexahedron, etc.)<br>
>>> PetscSection s;<br>
>>> KSP ksp;//Krylov Sub-Space (linear solver object)<br>
>>> Mat K,//Global stiffness matrix (Square, assume unsymmetric).<br>
>>> KE,//Element stiffness matrix (Square, assume unsymmetric).<br>
>>> matC;//Constitutive matrix.<br>
>>> Vec XYZ,//Coordinate vector, contains spatial locations of mesh's vertices (NOTE: This vector self-destroys!).<br>
>>> U,//Displacement vector.<br>
>>> F;//Load Vector.<br>
>>> PetscViewer XYZviewer,//Viewer object to output mesh to ASCII format.<br>
>>> XYZpUviewer; //Viewer object to output displacements to ASCII format.<br>
>>> PetscBool interpolate = PETSC_TRUE,//Instructs Gmsh importer whether to generate faces and edges (Needed when using P2 or higher elements).<br>
>>> useCone = PETSC_TRUE,//Instructs "DMPlexGetTransitiveClosure()" whether to extract the closure or the star.<br>
>>> dirichletBC = PETSC_FALSE,//For use when assembling the K matrix.<br>
>>> neumannBC = PETSC_FALSE,//For use when assembling the F vector.<br>
>>> saveASCII = PETSC_FALSE,//Whether to save results in ASCII format.<br>
>>> saveVTK = PETSC_FALSE;//Whether to save results as VTK format.<br>
>>> PetscInt nc,//number of cells. (PETSc lingo for "elements")<br>
>>> nv,//number of vertices. (PETSc lingo for "nodes")<br>
>>> nf,//number of faces. (PETSc lingo for "surfaces")<br>
>>> ne,//number of edges. (PETSc lingo for "lines")<br>
>>> pStart,//starting LI of global elements.<br>
>>> pEnd,//ending LI of all elements.<br>
>>> cStart,//starting LI for cells global arrangement.<br>
>>> cEnd,//ending LI for cells in global arrangement.<br>
>>> vStart,//starting LI for vertices in global arrangement.<br>
>>> vEnd,//ending LI for vertices in global arrangement.<br>
>>> fStart,//starting LI for faces in global arrangement.<br>
>>> fEnd,//ending LI for faces in global arrangement.<br>
>>> eStart,//starting LI for edges in global arrangement.<br>
>>> eEnd,//ending LI for edges in global arrangement.<br>
>>> sizeK,//Size of the element stiffness matrix.<br>
>>> ii,jj,kk,//Dedicated looping variables.<br>
>>> indexXYZ,//Variable to access the elements of XYZ vector.<br>
>>> indexK,//Variable to access the elements of the U and F vectors (can reference rows and colums of K matrix.)<br>
>>> *closure = PETSC_NULL,//Pointer to the closure elements of a cell.<br>
>>> size_closure,//Size of the closure of a cell.<br>
>>> dim,//Dimension of the mesh.<br>
>>> //*edof,//Linear indices of dof's inside the K matrix.<br>
>>> dof = 3,//Degrees of freedom per node.<br>
>>> cells=0, edges=0, vertices=0, faces=0,//Topology counters when looping through cells.<br>
>>> pinXcode=10, pinZcode=11,forceZcode=12;//Gmsh codes to extract relevant "Face Sets."<br>
>>> PetscReal //*x_el,//Pointer to a vector that will store the x-coordinates of an element's vertices.<br>
>>> //*y_el,//Pointer to a vector that will store the y-coordinates of an element's vertices.<br>
>>> //*z_el,//Pointer to a vector that will store the z-coordinates of an element's vertices.<br>
>>> *xyz_el,//Pointer to xyz array in the XYZ vector.<br>
>>> traction = -10,<br>
>>> *KEdata,<br>
>>> t1,t2; //time keepers.<br>
>>> const char *gmshfile = "TopOptmeshfine2.msh";//Name of gmsh file to import.<br>
>>><br>
>>> ierr = PetscInitialize(&argc,&args,NULL,help); if(ierr) return ierr; //And the machine shall work....<br>
>>><br>
>>> //MESH IMPORT=================================================================<br>
>>> //IMPORTANT NOTE: Gmsh only creates "cells" and "vertices," it does not create the "faces" or "edges."<br>
>>> //Gmsh probably can generate them, must figure out how to.<br>
>>> t1 = MPI_Wtime();<br>
>>> ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,gmshfile,interpolate,&dm); CHKERRQ(ierr);//Read Gmsh file and generate the DMPlex.<br>
>>> ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr);//1-D, 2-D, or 3-D<br>
>>> ierr = DMPlexGetChart(dm, &pStart, &pEnd); CHKERRQ(ierr);//Extracts linear indices of cells, vertices, faces, and edges.<br>
>>> ierr = DMGetCoordinatesLocal(dm,&XYZ); CHKERRQ(ierr);//Extracts coordinates from mesh.(Contiguous storage: [x0,y0,z0,x1,y1,z1,...])<br>
>>> ierr = VecGetArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's data.<br>
>>> t2 = MPI_Wtime();<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"Mesh Import time: %10f\n",t2-t1);<br>
>>> DMView(dm,PETSC_VIEWER_STDOUT_WORLD);<br>
>>><br>
>>> //IMPORTANT NOTE: PETSc assumes that vertex-cell meshes are 2D even if they were 3D, so its ordering changes.<br>
>>> //Cells remain at height 0, but vertices move to height 1 from height 3. To prevent this from becoming an issue<br>
>>> //the "interpolate" variable is set to PETSC_TRUE to tell the mesh importer to generate faces and edges.<br>
>>> //PETSc, therefore, technically does additional meshing. Gotta figure out how to get this from Gmsh directly.<br>
>>> ierr = DMPlexGetDepthStratum(dm,3, &cStart, &cEnd);//Get LI of cells.<br>
>>> ierr = DMPlexGetDepthStratum(dm,2, &fStart, &fEnd);//Get LI of faces<br>
>>> ierr = DMPlexGetDepthStratum(dm,1, &eStart, &eEnd);//Get LI of edges.<br>
>>> ierr = DMPlexGetDepthStratum(dm,0, &vStart, &vEnd);//Get LI of vertices.<br>
>>> ierr = DMGetStratumSize(dm,"depth", 3, &nc);//Get number of cells.<br>
>>> ierr = DMGetStratumSize(dm,"depth", 2, &nf);//Get number of faces.<br>
>>> ierr = DMGetStratumSize(dm,"depth", 1, &ne);//Get number of edges.<br>
>>> ierr = DMGetStratumSize(dm,"depth", 0, &nv);//Get number of vertices.<br>
>>> /*<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"global start = %10d\t global end = %10d\n",pStart,pEnd);<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"#cells = %10d\t i = %10d\t i < %10d\n",nc,cStart,cEnd);<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"#faces = %10d\t i = %10d\t i < %10d\n",nf,fStart,fEnd);<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"#edges = %10d\t i = %10d\t i < %10d\n",ne,eStart,eEnd);<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"#vertices = %10d\t i = %10d\t i < %10d\n",nv,vStart,vEnd);<br>
>>> */<br>
>>> //MESH IMPORT=================================================================<br>
>>><br>
>>> //NOTE: This section extremely hardcoded right now.<br>
>>> //Current setup would only support P1 meshes.<br>
>>> //MEMORY ALLOCATION ==========================================================<br>
>>> ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s); CHKERRQ(ierr);<br>
>>> //The chart is akin to a contiguous memory storage allocation. Each chart entry is associated<br>
>>> //with a "thing," could be a vertex, face, cell, or edge, or anything else.<br>
>>> ierr = PetscSectionSetChart(s, pStart, pEnd); CHKERRQ(ierr);<br>
>>> //For each "thing" in the chart, additional room can be made. This is helpful for associating<br>
>>> //nodes to multiple degrees of freedom. These commands help associate nodes with<br>
>>> for(ii = cStart; ii < cEnd; ii++){//Cell loop.<br>
>>> ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with cells.<br>
>>> for(ii = fStart; ii < fEnd; ii++){//Face loop.<br>
>>> ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with faces.<br>
>>> for(ii = vStart; ii < vEnd; ii++){//Vertex loop.<br>
>>> ierr = PetscSectionSetDof(s, ii, dof);CHKERRQ(ierr);}//Sets x, y, and z displacements as dofs.<br>
>>> for(ii = eStart; ii < eEnd; ii++){//Edge loop<br>
>>> ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with edges.<br>
>>> ierr = PetscSectionSetUp(s); CHKERRQ(ierr);<br>
>>> ierr = PetscSectionGetStorageSize(s,&sizeK);CHKERRQ(ierr);//Determine the size of the global stiffness matrix.<br>
>>> ierr = DMSetLocalSection(dm,s); CHKERRQ(ierr);//Associate the PetscSection with the DM object.<br>
>>> //PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)<br>
>>> //ierr = DMCreateGlobalVector(dm,&U); CHKERRQ(ierr);<br>
>>> PetscSectionDestroy(&s);<br>
>>> //PetscPrintf(PETSC_COMM_WORLD,"sizeK = %10d\n",sizeK);<br>
>>><br>
>>> //OBJECT SETUP================================================================<br>
>>> //Global stiffness matrix.<br>
>>> //PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)<br>
>>><br>
>>> //This makes the loop fast.<br>
>>> ierr = DMCreateMatrix(dm,&K);<br>
>>><br>
>>> //This makes the loop uber slow.<br>
>>> //ierr = MatCreate(PETSC_COMM_WORLD,&K); CHKERRQ(ierr);<br>
>>> //ierr = MatSetType(K,MATAIJ); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.<br>
>>> //ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,sizeK,sizeK); CHKERRQ(ierr);<br>
>>> //ierr = MatSetUp(K); CHKERRQ(ierr);<br>
>>><br>
>>> //Displacement vector.<br>
>>> ierr = VecCreate(PETSC_COMM_WORLD,&U); CHKERRQ(ierr);<br>
>>> ierr = VecSetType(U,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.<br>
>>> ierr = VecSetSizes(U,PETSC_DECIDE,sizeK); CHKERRQ(ierr);<br>
>>><br>
>>> //Load vector.<br>
>>> ierr = VecCreate(PETSC_COMM_WORLD,&F); CHKERRQ(ierr);<br>
>>> ierr = VecSetType(F,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.<br>
>>> ierr = VecSetSizes(F,PETSC_DECIDE,sizeK); CHKERRQ(ierr);<br>
>>> //OBJECT SETUP================================================================<br>
>>><br>
>>> //WARNING: This loop is currently hardcoded for P1 elements only! Must Figure<br>
>>> //out a clever way to modify to accomodate Pn (n>1) elements.<br>
>>><br>
>>> //BEGIN GLOBAL STIFFNESS MATRIX BUILDER=======================================<br>
>>> t1 = MPI_Wtime();<br>
>>><br>
>>> //PREALLOCATIONS==============================================================<br>
>>> ierr = ConstitutiveMatrix(&matC,"isotropic",0); CHKERRQ(ierr);<br>
>>> struct preKE preKEtetra4;<br>
>>> ierr = InitializeKEpreallocation(&preKEtetra4,"tetra4"); CHKERRQ(ierr);<br>
>>> ierr = MatCreate(PETSC_COMM_WORLD,&KE); CHKERRQ(ierr); //SEQUENTIAL<br>
>>> ierr = MatSetSizes(KE,PETSC_DECIDE,PETSC_DECIDE,12,12); CHKERRQ(ierr); //SEQUENTIAL<br>
>>> ierr = MatSetType(KE,MATDENSE); CHKERRQ(ierr); //SEQUENTIAL<br>
>>> ierr = MatSetUp(KE); CHKERRQ(ierr);<br>
>>> PetscReal x_tetra4[4], y_tetra4[4],z_tetra4[4],<br>
>>> x_hex8[8], y_hex8[8],z_hex8[8],<br>
>>> *x,*y,*z;<br>
>>> PetscInt *EDOF,edof_tetra4[12],edof_hex8[24];<br>
>>> DMPolytopeType previous = DM_POLYTOPE_UNKNOWN;<br>
>>> //PREALLOCATIONS==============================================================<br>
>>><br>
>>><br>
>>><br>
>>> for(ii=cStart;ii<cEnd;ii++){//loop through cells.<br>
>>> ierr = DMPlexGetTransitiveClosure(dm, ii, useCone, &size_closure, &closure); CHKERRQ(ierr);<br>
>>> ierr = DMPlexGetCellType(dm, ii, &celltype); CHKERRQ(ierr);<br>
>>> //IMPORTANT NOTE: MOST OF THIS LOOP SHOULD BE INCLUDED IN THE KE3D function.<br>
>>> if(previous != celltype){<br>
>>> //PetscPrintf(PETSC_COMM_WORLD,"run \n");<br>
>>> if(celltype == DM_POLYTOPE_TETRAHEDRON){<br>
>>> x = x_tetra4;<br>
>>> y = y_tetra4;<br>
>>> z = z_tetra4;<br>
>>> EDOF = edof_tetra4;<br>
>>> }//end if.<br>
>>> else if(celltype == DM_POLYTOPE_HEXAHEDRON){<br>
>>> x = x_hex8;<br>
>>> y = y_hex8;<br>
>>> z = z_hex8;<br>
>>> EDOF = edof_hex8;<br>
>>> }//end else if.<br>
>>> }<br>
>>> previous = celltype;<br>
>>><br>
>>> //PetscPrintf(PETSC_COMM_WORLD,"Cell # %4i\t",ii);<br>
>>> cells=0;<br>
>>> edges=0;<br>
>>> vertices=0;<br>
>>> faces=0;<br>
>>> kk = 0;<br>
>>> for(jj=0;jj<(2*size_closure);jj+=2){//Scan the closure of the current cell.<br>
>>> //Use information from the DM's strata to determine composition of cell_ii.<br>
>>> if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices.<br>
>>> //PetscPrintf(PETSC_COMM_WORLD,"%5i\t",closure[jj]);<br>
>>> indexXYZ = dim*(closure[jj]-vStart);//Linear index of x-coordinate in the xyz_el array.<br>
>>><br>
>>> *(x+vertices) = xyz_el[indexXYZ];<br>
>>> *(y+vertices) = xyz_el[indexXYZ+1];//Extract Y-coordinates of the current vertex.<br>
>>> *(z+vertices) = xyz_el[indexXYZ+2];//Extract Y-coordinates of the current vertex.<br>
>>> *(EDOF + kk) = indexXYZ;<br>
>>> *(EDOF + kk+1) = indexXYZ+1;<br>
>>> *(EDOF + kk+2) = indexXYZ+2;<br>
>>> kk+=3;<br>
>>> vertices++;//Update vertex counter.<br>
>>> }//end if<br>
>>> else if(eStart <= closure[jj] && closure[jj]< eEnd){//Check for edge ID's<br>
>>> edges++;<br>
>>> }//end else ifindexK<br>
>>> else if(fStart <= closure[jj] && closure[jj]< fEnd){//Check for face ID's<br>
>>> faces++;<br>
>>> }//end else if<br>
>>> else if(cStart <= closure[jj] && closure[jj]< cEnd){//Check for cell ID's<br>
>>> cells++;<br>
>>> }//end else if<br>
>>> }//end "jj" loop.<br>
>>> ierr = tetra4(x,y,z,&preKEtetra4,&matC,&KE); CHKERRQ(ierr); //Generate the element stiffness matrix for this cell.<br>
>>> ierr = MatDenseGetArray(KE,&KEdata); CHKERRQ(ierr);<br>
>>> ierr = MatSetValues(K,12,EDOF,12,EDOF,KEdata,ADD_VALUES); CHKERRQ(ierr);//WARNING: HARDCODED FOR TETRAHEDRAL P1 ELEMENTS ONLY !!!!!!!!!!!!!!!!!!!!!!!<br>
>>> ierr = MatDenseRestoreArray(KE,&KEdata); CHKERRQ(ierr);<br>
>>> ierr = DMPlexRestoreTransitiveClosure(dm, ii,useCone, &size_closure, &closure); CHKERRQ(ierr);<br>
>>> }//end "ii" loop.<br>
>>> ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> //ierr = MatView(K,PETSC_VIEWER_DRAW_WORLD); CHKERRQ(ierr);<br>
>>> //END GLOBAL STIFFNESS MATRIX BUILDER===========================================<br>
>>> t2 = MPI_Wtime();<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"K build time: %10f\n",t2-t1);<br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>><br>
>>> t1 = MPI_Wtime();<br>
>>> //BEGIN BOUNDARY CONDITION ENFORCEMENT==========================================<br>
>>> IS TrianglesIS, physicalsurfaceID;//, VerticesIS;<br>
>>> PetscInt numsurfvals,<br>
>>> //numRows,<br>
>>> dof_offset,numTri;<br>
>>> const PetscInt *surfvals,<br>
>>> //*pinZID,<br>
>>> *TriangleID;<br>
>>> PetscScalar diag =1;<br>
>>> PetscReal area,force;<br>
>>> //NOTE: Petsc can read/assign labels. Eeach label may posses multiple "values."<br>
>>> //These values act as tags within a tag.<br>
>>> //IMPORTANT NOTE: The below line needs a safety. If a mesh that does not feature<br>
>>> //face sets is imported, the code in its current state will crash!!!. This is currently<br>
>>> //hardcoded for the test mesh.<br>
>>> ierr = DMGetLabel(dm, "Face Sets", &physicalgroups); CHKERRQ(ierr);//Inspects Physical surface groups defined by gmsh (if any).<br>
>>> ierr = DMLabelGetValueIS(physicalgroups, &physicalsurfaceID); CHKERRQ(ierr);//Gets the physical surface ID's defined in gmsh (as specified in the .geo file).<br>
>>> ierr = ISGetIndices(physicalsurfaceID,&surfvals); CHKERRQ(ierr);//Get a pointer to the actual surface values.<br>
>>> ierr = DMLabelGetNumValues(physicalgroups, &numsurfvals); CHKERRQ(ierr);//Gets the number of different values that the label assigns.<br>
>>> for(ii=0;ii<numsurfvals;ii++){//loop through the values under the label.<br>
>>> //PetscPrintf(PETSC_COMM_WORLD,"Values = %5i\n",surfvals[ii]);<br>
>>> //PROBLEM: The surface values are hardcoded in the gmsh file. We need to adopt standard "codes"<br>
>>> //that we can give to users when they make their meshes so that this code recognizes the Type<br>
>>> // of boundary conditions that are to be imposed.<br>
>>> if(surfvals[ii] == pinXcode){<br>
>>> dof_offset = 0;<br>
>>> dirichletBC = PETSC_TRUE;<br>
>>> }//end if.<br>
>>> else if(surfvals[ii] == pinZcode){<br>
>>> dof_offset = 2;<br>
>>> dirichletBC = PETSC_TRUE;<br>
>>> }//end else if.<br>
>>> else if(surfvals[ii] == forceZcode){<br>
>>> dof_offset = 2;<br>
>>> neumannBC = PETSC_TRUE;<br>
>>> }//end else if.<br>
>>><br>
>>> ierr = DMLabelGetStratumIS(physicalgroups, surfvals[ii], &TrianglesIS); CHKERRQ(ierr);//Get the ID's (as an IS) of the surfaces belonging to value 11.<br>
>>> //PROBLEM: DMPlexGetConeRecursiveVertices returns an array with repeated node ID's. For each repetition, the lines that enforce BC's unnecessarily re-run.<br>
>>> ierr = ISGetSize(TrianglesIS,&numTri); CHKERRQ(ierr);<br>
>>> ierr = ISGetIndices(TrianglesIS,&TriangleID); CHKERRQ(ierr);//Get a pointer to the actual surface values.<br>
>>> for(kk=0;kk<numTri;kk++){<br>
>>> ierr = DMPlexGetTransitiveClosure(dm, TriangleID[kk], useCone, &size_closure, &closure); CHKERRQ(ierr);<br>
>>> if(neumannBC){<br>
>>> ierr = DMPlexComputeCellGeometryFVM(dm, TriangleID[kk], &area,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);<br>
>>> force = traction*area/3;//WARNING: The 3 here is hardcoded for a purely tetrahedral mesh only!!!!!!!!!!<br>
>>> }<br>
>>> for(jj=0;jj<(2*size_closure);jj+=2){<br>
>>> //PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])<br>
>>> if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices.<br>
>>> indexK = dof*(closure[jj] - vStart) + dof_offset; //Compute the dof ID's in the K matrix.<br>
>>> if(dirichletBC){//Boundary conditions requiring an edit of K matrix.<br>
>>> ierr = MatZeroRows(K,1,&indexK,diag,NULL,NULL); CHKERRQ(ierr);<br>
>>> }//end if.<br>
>>> else if(neumannBC){//Boundary conditions requiring an edit of RHS vector.<br>
>>> ierr = VecSetValue(F,indexK,force,ADD_VALUES); CHKERRQ(ierr);<br>
>>> }// end else if.<br>
>>> }//end if.<br>
>>> }//end "jj" loop.<br>
>>> ierr = DMPlexRestoreTransitiveClosure(dm, closure[jj],useCone, &size_closure, &closure); CHKERRQ(ierr);<br>
>>> }//end "kk" loop.<br>
>>> ierr = ISRestoreIndices(TrianglesIS,&TriangleID); CHKERRQ(ierr);<br>
>>><br>
>>> /*<br>
>>> 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.<br>
>>> ierr = ISGetSize(VerticesIS,&numRows); CHKERRQ(ierr);//Get number of flagged vertices (this includes repeated indices for faces that share nodes).<br>
>>> ierr = ISGetIndices(VerticesIS,&pinZID); CHKERRQ(ierr);//Get a pointer to the actual surface values.<br>
>>> if(dirichletBC){//Boundary conditions requiring an edit of K matrix.<br>
>>> for(kk=0;kk<numRows;kk++){<br>
>>> 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.)<br>
>>> ierr = MatZeroRows(K,1,&indexK,diag,NULL,NULL); CHKERRQ(ierr);<br>
>>> }//end "kk" loop.<br>
>>> }//end if.<br>
>>> else if(neumannBC){//Boundary conditions requiring an edit of RHS vector.<br>
>>> for(kk=0;kk<numRows;kk++){<br>
>>> indexK = 3*(pinZID[kk] - vStart) + dof_offset;<br>
>>> ierr = VecSetValue(F,indexK,traction,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> }//end "kk" loop.<br>
>>> }// end else if.<br>
>>> ierr = ISRestoreIndices(VerticesIS,&pinZID); CHKERRQ(ierr);<br>
>>> */<br>
>>> dirichletBC = PETSC_FALSE;<br>
>>> neumannBC = PETSC_FALSE;<br>
>>> }//end "ii" loop.<br>
>>> ierr = ISRestoreIndices(physicalsurfaceID,&surfvals); CHKERRQ(ierr);<br>
>>> //ierr = ISRestoreIndices(VerticesIS,&pinZID); CHKERRQ(ierr);<br>
>>> ierr = ISDestroy(&physicalsurfaceID); CHKERRQ(ierr);<br>
>>> //ierr = ISDestroy(&VerticesIS); CHKERRQ(ierr);<br>
>>> ierr = ISDestroy(&TrianglesIS); CHKERRQ(ierr);<br>
>>> //END BOUNDARY CONDITION ENFORCEMENT============================================<br>
>>> t2 = MPI_Wtime();<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"BC imposition time: %10f\n",t2-t1);<br>
>>><br>
>>> /*<br>
>>> PetscInt kk = 0;<br>
>>> for(ii=vStart;ii<vEnd;ii++){<br>
>>> kk++;<br>
>>> 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]);<br>
>>> }// end "ii" loop.<br>
>>> */<br>
>>><br>
>>> t1 = MPI_Wtime();<br>
>>> //SOLVER========================================================================<br>
>>> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);<br>
>>> ierr = KSPSetOperators(ksp,K,K); CHKERRQ(ierr);<br>
>>> ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);<br>
>>> ierr = KSPSolve(ksp,F,U); CHKERRQ(ierr);<br>
>>> t2 = MPI_Wtime();<br>
>>> //ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br>
>>> //SOLVER========================================================================<br>
>>> t2 = MPI_Wtime();<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"Solver time: %10f\n",t2-t1);<br>
>>> ierr = VecRestoreArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's data.<br>
>>><br>
>>> //BEGIN MAX/MIN DISPLACEMENTS===================================================<br>
>>> IS ISux,ISuy,ISuz;<br>
>>> Vec UX,UY,UZ;<br>
>>> PetscReal UXmax,UYmax,UZmax,UXmin,UYmin,UZmin;<br>
>>> ierr = ISCreateStride(PETSC_COMM_WORLD,nv,0,3,&ISux); CHKERRQ(ierr);<br>
>>> ierr = ISCreateStride(PETSC_COMM_WORLD,nv,1,3,&ISuy); CHKERRQ(ierr);<br>
>>> ierr = ISCreateStride(PETSC_COMM_WORLD,nv,2,3,&ISuz); CHKERRQ(ierr);<br>
>>><br>
>>> //PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)<br>
>>> ierr = VecGetSubVector(U,ISux,&UX); CHKERRQ(ierr);<br>
>>> ierr = VecGetSubVector(U,ISuy,&UY); CHKERRQ(ierr);<br>
>>> ierr = VecGetSubVector(U,ISuz,&UZ); CHKERRQ(ierr);<br>
>>><br>
>>> //PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)<br>
>>> ierr = VecMax(UX,PETSC_NULL,&UXmax); CHKERRQ(ierr);<br>
>>> ierr = VecMax(UY,PETSC_NULL,&UYmax); CHKERRQ(ierr);<br>
>>> ierr = VecMax(UZ,PETSC_NULL,&UZmax); CHKERRQ(ierr);<br>
>>><br>
>>> ierr = VecMin(UX,PETSC_NULL,&UXmin); CHKERRQ(ierr);<br>
>>> ierr = VecMin(UY,PETSC_NULL,&UYmin); CHKERRQ(ierr);<br>
>>> ierr = VecMin(UZ,PETSC_NULL,&UZmin); CHKERRQ(ierr);<br>
>>><br>
>>> PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= ux <= %10f\n",UXmin,UXmax);<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= uy <= %10f\n",UYmin,UYmax);<br>
>>> PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= uz <= %10f\n",UZmin,UZmax);<br>
>>><br>
>>><br>
>>><br>
>>><br>
>>> //BEGIN OUTPUT SOLUTION=========================================================<br>
>>> if(saveASCII){<br>
>>> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"XYZ.txt",&XYZviewer);<br>
>>> ierr = VecView(XYZ,XYZviewer); CHKERRQ(ierr);<br>
>>> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"U.txt",&XYZpUviewer);<br>
>>> ierr = VecView(U,XYZpUviewer); CHKERRQ(ierr);<br>
>>> PetscViewerDestroy(&XYZviewer); PetscViewerDestroy(&XYZpUviewer);<br>
>>><br>
>>> }//end if.<br>
>>> if(saveVTK){<br>
>>> const char *meshfile = "starting_mesh.vtk",<br>
>>> *deformedfile = "deformed_mesh.vtk";<br>
>>> ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,meshfile,FILE_MODE_WRITE,&XYZviewer); CHKERRQ(ierr);<br>
>>> //PetscErrorCode DMSetAuxiliaryVec(DM dm, DMLabel label, PetscInt value, Vec aux)<br>
>>> DMLabel UXlabel,UYlabel, UZlabel;<br>
>>> //PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)<br>
>>> ierr = DMLabelCreate(PETSC_COMM_WORLD, "X-Displacement", &UXlabel); CHKERRQ(ierr);<br>
>>> ierr = DMLabelCreate(PETSC_COMM_WORLD, "Y-Displacement", &UYlabel); CHKERRQ(ierr);<br>
>>> ierr = DMLabelCreate(PETSC_COMM_WORLD, "Z-Displacement", &UZlabel); CHKERRQ(ierr);<br>
>>> ierr = DMSetAuxiliaryVec(dm,UXlabel, 1, UX); CHKERRQ(ierr);<br>
>>> ierr = DMSetAuxiliaryVec(dm,UYlabel, 1, UY); CHKERRQ(ierr);<br>
>>> ierr = DMSetAuxiliaryVec(dm,UZlabel, 1, UZ); CHKERRQ(ierr);<br>
>>> //PetscErrorCode PetscViewerVTKAddField(PetscViewer viewer,PetscObject dm,PetscErrorCode (*PetscViewerVTKWriteFunction)(PetscObject,PetscViewer),PetscInt fieldnum,PetscViewerVTKFieldType fieldtype,PetscBool checkdm,PetscObject vec)<br>
>>><br>
>>><br>
>>><br>
>>> //ierr = PetscViewerVTKAddField(XYZviewer, dm,PetscErrorCode (*PetscViewerVTKWriteFunction)(Vec,PetscViewer),PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,UX);<br>
>>> ierr = PetscViewerVTKAddField(XYZviewer, (PetscObject)dm,&PetscViewerVTKWriteFunction,PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,(PetscObject)UX);<br>
>>><br>
>>><br>
>>> ierr = DMPlexVTKWriteAll((PetscObject)dm, XYZviewer); CHKERRQ(ierr);<br>
>>> ierr = VecAXPY(XYZ,1,U); CHKERRQ(ierr);//Add displacement field to the mesh coordinates to deform.<br>
>>> ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,deformedfile,FILE_MODE_WRITE,&XYZpUviewer); CHKERRQ(ierr);<br>
>>> ierr = DMPlexVTKWriteAll((PetscObject)dm, XYZpUviewer); CHKERRQ(ierr);//<br>
>>> PetscViewerDestroy(&XYZviewer); PetscViewerDestroy(&XYZpUviewer);<br>
>>><br>
>>> }//end else if.<br>
>>> else{<br>
>>> ierr = PetscPrintf(PETSC_COMM_WORLD,"No output format specified! Files not saved.\n"); CHKERRQ(ierr);<br>
>>> }//end else.<br>
>>><br>
>>><br>
>>> //END OUTPUT SOLUTION===========================================================<br>
>>> VecDestroy(&UX); ISDestroy(&ISux);<br>
>>> VecDestroy(&UY); ISDestroy(&ISuy);<br>
>>> VecDestroy(&UZ); ISDestroy(&ISuz);<br>
>>> //END MAX/MIN DISPLACEMENTS=====================================================<br>
>>><br>
>>> //CLEANUP=====================================================================<br>
>>> DMDestroy(&dm);<br>
>>> KSPDestroy(&ksp);<br>
>>> MatDestroy(&K); MatDestroy(&KE); MatDestroy(&matC); //MatDestroy(preKEtetra4.matB); MatDestroy(preKEtetra4.matBTCB);<br>
>>> VecDestroy(&U); VecDestroy(&F);<br>
>>><br>
>>> //DMLabelDestroy(&physicalgroups);//Destroyig the DM destroys the label.<br>
>>> //CLEANUP=====================================================================<br>
>>> //PetscErrorCode PetscMallocDump(FILE *fp)<br>
>>> //ierr = PetscMallocDump(NULL);<br>
>>> return PetscFinalize();//And the machine shall rest....<br>
>>> }//end main.<br>
>>><br>
>>> PetscErrorCode tetra4(PetscScalar* X,PetscScalar* Y, PetscScalar* Z,struct preKE *P, Mat* matC, Mat* KE){<br>
>>> //INPUTS:<br>
>>> //X: Global X coordinates of the elemental nodes.<br>
>>> //Y: Global Y coordinates of the elemental nodes.<br>
>>> //Z: Global Z coordinates of the elemental nodes.<br>
>>> //J: Jacobian matrix.<br>
>>> //invJ: Inverse Jacobian matrix.<br>
>>> PetscErrorCode ierr;<br>
>>> //For current quadrature point, get dPsi/dXi_i Xi_i = {Xi,Eta,Zeta}<br>
>>> /*<br>
>>> P->dPsidXi[0] = +1.; P->dPsidEta[0] = 0.0; P->dPsidZeta[0] = 0.0;<br>
>>> P->dPsidXi[1] = 0.0; P->dPsidEta[1] = +1.; P->dPsidZeta[1] = 0.0;<br>
>>> P->dPsidXi[2] = 0.0; P->dPsidEta[2] = 0.0; P->dPsidZeta[2] = +1.;<br>
>>> P->dPsidXi[3] = -1.; P->dPsidEta[3] = -1.; P->dPsidZeta[3] = -1.;<br>
>>> */<br>
>>> //Populate the Jacobian matrix.<br>
>>> P->J[0][0] = X[0] - X[3];<br>
>>> P->J[0][1] = Y[0] - Y[3];<br>
>>> P->J[0][2] = Z[0] - Z[3];<br>
>>> P->J[1][0] = X[1] - X[3];<br>
>>> P->J[1][1] = Y[1] - Y[3];<br>
>>> P->J[1][2] = Z[1] - Z[3];<br>
>>> P->J[2][0] = X[2] - X[3];<br>
>>> P->J[2][1] = Y[2] - Y[3];<br>
>>> P->J[2][2] = Z[2] - Z[3];<br>
>>><br>
>>> //Determinant of the 3x3 Jacobian. (Expansion along 1st row).<br>
>>> P->minor00 = P->J[1][1]*P->J[2][2] - P->J[2][1]*P->J[1][2];//Reuse when finding InvJ.<br>
>>> P->minor01 = P->J[1][0]*P->J[2][2] - P->J[2][0]*P->J[1][2];//Reuse when finding InvJ.<br>
>>> P->minor02 = P->J[1][0]*P->J[2][1] - P->J[2][0]*P->J[1][1];//Reuse when finding InvJ.<br>
>>> P->detJ = P->J[0][0]*P->minor00 - P->J[0][1]*P->minor01 + P->J[0][2]*P->minor02;<br>
>>> //Inverse of the 3x3 Jacobian<br>
>>> P->invJ[0][0] = +P->minor00/P->detJ;//Reuse precomputed minor.<br>
>>> P->invJ[0][1] = -(P->J[0][1]*P->J[2][2] - P->J[0][2]*P->J[2][1])/P->detJ;<br>
>>> P->invJ[0][2] = +(P->J[0][1]*P->J[1][2] - P->J[1][1]*P->J[0][2])/P->detJ;<br>
>>> P->invJ[1][0] = -P->minor01/P->detJ;//Reuse precomputed minor.<br>
>>> P->invJ[1][1] = +(P->J[0][0]*P->J[2][2] - P->J[0][2]*P->J[2][0])/P->detJ;<br>
>>> P->invJ[1][2] = -(P->J[0][0]*P->J[1][2] - P->J[1][0]*P->J[0][2])/P->detJ;<br>
>>> P->invJ[2][0] = +P->minor02/P->detJ;//Reuse precomputed minor.<br>
>>> P->invJ[2][1] = -(P->J[0][0]*P->J[2][1] - P->J[0][1]*P->J[2][0])/P->detJ;<br>
>>> P->invJ[2][2] = +(P->J[0][0]*P->J[1][1] - P->J[0][1]*P->J[1][0])/P->detJ;<br>
>>><br>
>>> //*****************STRAIN MATRIX (B)**************************************<br>
>>> for(P->m=0;P->m<P->N;P->m++){//Scan all shape functions.<br>
>>><br>
>>> P->x_in = 0 + P->m*3;//Every 3rd column starting at 0<br>
>>> P->y_in = P->x_in +1;//Every 3rd column starting at 1<br>
>>> P->z_in = P->y_in +1;//Every 3rd column starting at 2<br>
>>><br>
>>> 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];<br>
>>> 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];<br>
>>> 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];<br>
>>><br>
>>> P->dX[1] = P->dZ[0]; P->dX[2] = P->dY[0];<br>
>>> P->dY[1] = P->dZ[0]; P->dY[2] = P->dX[0];<br>
>>> P->dZ[1] = P->dX[0]; P->dZ[2] = P->dY[0];<br>
>>><br>
>>> ierr = MatSetValues(P->matB,3,P->x_insert,1,&(P->x_in),P->dX,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> ierr = MatSetValues(P->matB,3,P->y_insert,1,&(P->y_in),P->dY,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> ierr = MatSetValues(P->matB,3,P->z_insert,1,&(P->z_in),P->dZ,INSERT_VALUES); CHKERRQ(ierr);<br>
>>><br>
>>> }//end "m" loop.<br>
>>> ierr = MatAssemblyBegin(P->matB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> ierr = MatAssemblyEnd(P->matB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> //*****************STRAIN MATRIX (B)**************************************<br>
>>><br>
>>> //Compute the matrix product B^t*C*B, scale it by the quadrature weights and add to KE.<br>
>>> P->weight = -P->detJ/6;<br>
>>><br>
>>> ierr = MatZeroEntries(*KE); CHKERRQ(ierr);<br>
>>> ierr = MatPtAP(*matC,P->matB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(P->matBTCB));CHKERRQ(ierr);<br>
>>> ierr = MatScale(P->matBTCB,P->weight); CHKERRQ(ierr);<br>
>>> ierr = MatAssemblyBegin(P->matBTCB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> ierr = MatAssemblyEnd(P->matBTCB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> ierr = MatAXPY(*KE,1,P->matBTCB,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);//Add contribution of current quadrature point to KE.<br>
>>><br>
>>> //ierr = MatPtAP(*matC,P->matB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,KE);CHKERRQ(ierr);<br>
>>> //ierr = MatScale(*KE,P->weight); CHKERRQ(ierr);<br>
>>><br>
>>> ierr = MatAssemblyBegin(*KE,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> ierr = MatAssemblyEnd(*KE,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>><br>
>>> //Cleanup<br>
>>> return ierr;<br>
>>> }//end tetra4.<br>
>>><br>
>>> PetscErrorCode ConstitutiveMatrix(Mat *matC,const char* type,PetscInt materialID){<br>
>>> PetscErrorCode ierr;<br>
>>> PetscBool isotropic = PETSC_FALSE,<br>
>>> orthotropic = PETSC_FALSE;<br>
>>> //PetscErrorCode PetscStrcmp(const char a[],const char b[],PetscBool *flg)<br>
>>> ierr = PetscStrcmp(type,"isotropic",&isotropic);<br>
>>> ierr = PetscStrcmp(type,"orthotropic",&orthotropic);<br>
>>> ierr = MatCreate(PETSC_COMM_WORLD,matC); CHKERRQ(ierr);<br>
>>> ierr = MatSetSizes(*matC,PETSC_DECIDE,PETSC_DECIDE,6,6); CHKERRQ(ierr);<br>
>>> ierr = MatSetType(*matC,MATAIJ); CHKERRQ(ierr);<br>
>>> ierr = MatSetUp(*matC); CHKERRQ(ierr);<br>
>>><br>
>>> if(isotropic){<br>
>>> PetscReal E,nu, M,L,vals[3];<br>
>>> switch(materialID){<br>
>>> case 0://Hardcoded properties for isotropic material #0<br>
>>> E = 200;<br>
>>> nu = 1./3;<br>
>>> break;<br>
>>> case 1://Hardcoded properties for isotropic material #1<br>
>>> E = 96;<br>
>>> nu = 1./3;<br>
>>> break;<br>
>>> }//end switch.<br>
>>> M = E/(2*(1+nu)),//Lame's constant 1 ("mu").<br>
>>> L = E*nu/((1+nu)*(1-2*nu));//Lame's constant 2 ("lambda").<br>
>>> //PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)<br>
>>> PetscInt idxn[3] = {0,1,2};<br>
>>> vals[0] = L+2*M; vals[1] = L; vals[2] = vals[1];<br>
>>> ierr = MatSetValues(*matC,1,&idxn[0],3,idxn,vals,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> vals[1] = vals[0]; vals[0] = vals[2];<br>
>>> ierr = MatSetValues(*matC,1,&idxn[1],3,idxn,vals,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> vals[2] = vals[1]; vals[1] = vals[0];<br>
>>> ierr = MatSetValues(*matC,1,&idxn[2],3,idxn,vals,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> ierr = MatSetValue(*matC,3,3,M,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> ierr = MatSetValue(*matC,4,4,M,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> ierr = MatSetValue(*matC,5,5,M,INSERT_VALUES); CHKERRQ(ierr);<br>
>>> }//end if.<br>
>>> /*<br>
>>> else if(orthotropic){<br>
>>> switch(materialID){<br>
>>> case 0:<br>
>>> break;<br>
>>> case 1:<br>
>>> break;<br>
>>> }//end switch.<br>
>>> }//end else if.<br>
>>> */<br>
>>> ierr = MatAssemblyBegin(*matC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> ierr = MatAssemblyEnd(*matC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
>>> //MatView(*matC,0);<br>
>>> return ierr;<br>
>>> }//End ConstitutiveMatrix<br>
>>><br>
>>> PetscErrorCode InitializeKEpreallocation(struct preKE *P,const char* type){<br>
>>> PetscErrorCode ierr;<br>
>>> PetscBool istetra4 = PETSC_FALSE,<br>
>>> ishex8 = PETSC_FALSE;<br>
>>> ierr = PetscStrcmp(type,"tetra4",&istetra4); CHKERRQ(ierr);<br>
>>> ierr = PetscStrcmp(type,"hex8",&ishex8); CHKERRQ(ierr);<br>
>>> if(istetra4){<br>
>>> P->sizeKE = 12;<br>
>>> P->N = 4;<br>
>>> }//end if.<br>
>>> else if(ishex8){<br>
>>> P->sizeKE = 24;<br>
>>> P->N = 8;<br>
>>> }//end else if.<br>
>>><br>
>>><br>
>>> P->x_insert[0] = 0; P->x_insert[1] = 3; P->x_insert[2] = 5;<br>
>>> P->y_insert[0] = 1; P->y_insert[1] = 4; P->y_insert[2] = 5;<br>
>>> P->z_insert[0] = 2; P->z_insert[1] = 3; P->z_insert[2] = 4;<br>
>>> //Allocate memory for the differentiated shape function vectors.<br>
>>> ierr = PetscMalloc1(P->N,&(P->dPsidXi)); CHKERRQ(ierr);<br>
>>> ierr = PetscMalloc1(P->N,&(P->dPsidEta)); CHKERRQ(ierr);<br>
>>> ierr = PetscMalloc1(P->N,&(P->dPsidZeta)); CHKERRQ(ierr);<br>
>>><br>
>>> P->dPsidXi[0] = +1.; P->dPsidEta[0] = 0.0; P->dPsidZeta[0] = 0.0;<br>
>>> P->dPsidXi[1] = 0.0; P->dPsidEta[1] = +1.; P->dPsidZeta[1] = 0.0;<br>
>>> P->dPsidXi[2] = 0.0; P->dPsidEta[2] = 0.0; P->dPsidZeta[2] = +1.;<br>
>>> P->dPsidXi[3] = -1.; P->dPsidEta[3] = -1.; P->dPsidZeta[3] = -1.;<br>
>>><br>
>>><br>
>>> //Strain matrix.<br>
>>> ierr = MatCreate(PETSC_COMM_WORLD,&(P->matB)); CHKERRQ(ierr);<br>
>>> ierr = MatSetSizes(P->matB,PETSC_DECIDE,PETSC_DECIDE,6,P->sizeKE); CHKERRQ(ierr);//Hardcoded<br>
>>> ierr = MatSetType(P->matB,MATAIJ); CHKERRQ(ierr);<br>
>>> ierr = MatSetUp(P->matB); CHKERRQ(ierr);<br>
>>><br>
>>> //Contribution matrix.<br>
>>> ierr = MatCreate(PETSC_COMM_WORLD,&(P->matBTCB)); CHKERRQ(ierr);<br>
>>> ierr = MatSetSizes(P->matBTCB,PETSC_DECIDE,PETSC_DECIDE,P->sizeKE,P->sizeKE); CHKERRQ(ierr);<br>
>>> ierr = MatSetType(P->matBTCB,MATAIJ); CHKERRQ(ierr);<br>
>>> ierr = MatSetUp(P->matBTCB); CHKERRQ(ierr);<br>
>>><br>
>>> //Element stiffness matrix.<br>
>>> //ierr = MatCreateSeqDense(PETSC_COMM_SELF,12,12,NULL,&KE); CHKERRQ(ierr); //PARALLEL<br>
>>><br>
>>> return ierr;<br>
>>> }<br>
</div>
</span></font></div>
</div>
</body>
</html>