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