Actual source code: cartesian.c

petsc-3.3-p7 2013-05-11
  1: #include <petsc-private/meshimpl.h>   /*I      "petscdmmesh.h"   I*/
  2: #include <CartesianSieve.hh>

  6: /*@C
  7:   DMCartesianGetMesh - Gets the internal mesh object

  9:   Not collective

 11:   Input Parameter:
 12: . dm - the mesh object

 14:   Output Parameter:
 15: . m - the internal mesh object

 17:   Level: advanced

 19: .seealso MeshCreate(), MeshCartesianSetMesh()
 20: @*/
 21: PetscErrorCode DMCartesianGetMesh(DM dm, ALE::Obj<ALE::CartesianMesh>& m)
 22: {
 23:   DM_Cartesian *c = (DM_Cartesian *) dm->data;

 27:   m = c->m;
 28:   return(0);
 29: }

 33: /*@C
 34:   DMCartesianSetMesh - Sets the internal mesh object

 36:   Not collective

 38:   Input Parameters:
 39: + mesh - the mesh object
 40: - m - the internal mesh object

 42:   Level: advanced

 44: .seealso MeshCreate(), MeshCartesianGetMesh()
 45: @*/
 46: PetscErrorCode DMCartesianSetMesh(DM dm, const ALE::Obj<ALE::CartesianMesh>& m)
 47: {
 48:   DM_Cartesian *c = (DM_Cartesian *) dm->data;

 52:   c->m = m;
 53:   return(0);
 54: }

 58: PetscErrorCode  DMDestroy_Cartesian(DM dm)
 59: {
 60:   DM_Cartesian *c = (DM_Cartesian *) dm->data;

 63:   c->m = PETSC_NULL;
 64:   return(0);
 65: }

 69: PetscErrorCode DMView_Cartesian_Ascii(const ALE::Obj<ALE::CartesianMesh>& mesh, PetscViewer viewer)
 70: {
 71:   PetscViewerFormat format;
 72:   PetscErrorCode    ierr;

 75:   PetscViewerGetFormat(viewer, &format);
 76:   if (format == PETSC_VIEWER_ASCII_VTK) {
 77: #if 0
 78:     VTKViewer::writeHeader(viewer);
 79:     VTKViewer::writeVertices(mesh, viewer);
 80:     VTKViewer::writeElements(mesh, viewer);
 81: #endif
 82:   } else {
 83:     int dim = mesh->getDimension();

 85:     PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
 86:     /* FIX: Need to globalize */
 87:     PetscViewerASCIIPrintf(viewer, "  %d vertices\n", mesh->getSieve()->getNumVertices());
 88:     PetscViewerASCIIPrintf(viewer, "  %d cells\n",    mesh->getSieve()->getNumCells());
 89:   }
 90:   PetscViewerFlush(viewer);
 91:   return(0);
 92: }

 94: PetscErrorCode DMView_Cartesian(DM dm, PetscViewer viewer)
 95: {
 96:   ALE::Obj<ALE::CartesianMesh> m;
 97:   PetscBool      iascii, isbinary, isdraw;

101:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
102:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
103:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);

105:   DMCartesianGetMesh(dm, m);
106:   if (iascii){
107:     DMView_Cartesian_Ascii(m, viewer);
108:   } else if (isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Binary viewer not implemented for Cartesian Mesh");
109:   else if (isdraw) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Draw viewer not implemented for Cartesian Mesh");
110:   else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
111:   return(0);
112: }

116: PetscErrorCode DMCreateInterpolation_Cartesian(DM fineMesh, DM coarseMesh, Mat *interpolation, Vec *scaling)
117: {
118:   ALE::Obj<ALE::CartesianMesh> coarse;
119:   ALE::Obj<ALE::CartesianMesh> fine;
120:   Mat                          P;
121:   PetscErrorCode               ierr;

124:   DMCartesianGetMesh(fineMesh,   fine);
125:   DMCartesianGetMesh(coarseMesh, coarse);
126: #if 0
127:   const ALE::Obj<ALE::Mesh::real_section_type>& coarseCoordinates = coarse->getRealSection("coordinates");
128:   const ALE::Obj<ALE::Mesh::real_section_type>& fineCoordinates   = fine->getRealSection("coordinates");
129:   const ALE::Obj<ALE::Mesh::label_sequence>&    vertices          = fine->depthStratum(0);
130:   const ALE::Obj<ALE::Mesh::real_section_type>& sCoarse           = coarse->getRealSection("default");
131:   const ALE::Obj<ALE::Mesh::real_section_type>& sFine             = fine->getRealSection("default");
132:   const ALE::Obj<ALE::Mesh::order_type>&        coarseOrder = coarse->getFactory()->getGlobalOrder(coarse, "default", sCoarse);
133:   const ALE::Obj<ALE::Mesh::order_type>&        fineOrder   = fine->getFactory()->getGlobalOrder(fine, "default", sFine);
134:   const int dim = coarse->getDimension();
135:   double *v0, *J, *invJ, detJ, *refCoords, *values;
136: #endif

138:   MatCreate(fine->comm(), &P);
139: #if 0
140:   MatSetSizes(P, sFine->size(), sCoarse->size(), PETSC_DETERMINE, PETSC_DETERMINE);
141:   MatSetFromOptions(P);
142:   PetscMalloc5(dim,double,&v0,dim*dim,double,&J,dim*dim,double,&invJ,dim,double,&refCoords,dim+1,double,&values);
143:   for(ALE::Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
144:     const ALE::Mesh::real_section_type::value_type *coords     = fineCoordinates->restrictPoint(*v_iter);
145:     const ALE::Mesh::point_type                     coarseCell = coarse->locatePoint(coords);

147:     coarse->computeElementGeometry(coarseCoordinates, coarseCell, v0, J, invJ, detJ);
148:     for(int d = 0; d < dim; ++d) {
149:       refCoords[d] = 0.0;
150:       for(int e = 0; e < dim; ++e) {
151:         refCoords[d] += invJ[d*dim+e]*(coords[e] - v0[e]);
152:       }
153:       refCoords[d] -= 1.0;
154:     }
155:     values[0] = 1.0/3.0 - (refCoords[0] + refCoords[1])/3.0;
156:     values[1] = 0.5*(refCoords[0] + 1.0);
157:     values[2] = 0.5*(refCoords[1] + 1.0);
158:     updateOperatorGeneral(P, fine, sFine, fineOrder, *v_iter, sCoarse, coarseOrder, coarseCell, values, INSERT_VALUES);
159:   }
160:   PetscFree5(v0,J,invJ,refCoords,values);
161: #endif
162:   MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
163:   MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
164:   *interpolation = P;
165:   return(0);
166: }

170: PetscErrorCode DMRefine_Cartesian(DM mesh, MPI_Comm comm, DM *refinedMesh)
171: {
172:   ALE::Obj<ALE::CartesianMesh> oldMesh;
173:   PetscErrorCode               ierr;

176:   DMCartesianGetMesh(mesh, oldMesh);
177:   DMCartesianCreate(comm, refinedMesh);
178: #if 0
179:   ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
180:   MeshCartesianSetMesh(*refinedMesh, newMesh);
181:   const ALE::Obj<ALE::CartesianMesh::real_section_type>& s = newMesh->getRealSection("default");

183:   newMesh->setDiscretization(oldMesh->getDiscretization());
184:   newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
185:   newMesh->setupField(s);
186: #else
187:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Not yet implemented");
188: #endif
189:   return(0);
190: }

194: PetscErrorCode DMCoarsen_Cartesian(DM mesh, MPI_Comm comm, DM *coarseMesh)
195: {
196:   ALE::Obj<ALE::CartesianMesh> oldMesh;
197:   PetscErrorCode               ierr;

200:   if (comm == MPI_COMM_NULL) comm = ((PetscObject)mesh)->comm;
201:   DMCartesianGetMesh(mesh, oldMesh);
202:   DMCartesianCreate(comm, coarseMesh);
203: #if 0
204:   ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
205:   MeshCartesianSetMesh(*coarseMesh, newMesh);
206:   const ALE::Obj<ALE::CartesianMesh::real_section_type>& s = newMesh->getRealSection("default");

208:   newMesh->setDiscretization(oldMesh->getDiscretization());
209:   newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
210:   newMesh->setupField(s);
211: #else
212:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Not yet implemented");
213: #endif
214:   return(0);
215: }

219: PetscErrorCode DMCartesianGetSectionReal(DM dm, const char name[], SectionReal *section)
220: {
221:   ALE::Obj<ALE::CartesianMesh> m;
222:   PetscErrorCode      ierr;

225:   DMCartesianGetMesh(dm, m);
226:   SectionRealCreate(m->comm(), section);
227:   PetscObjectSetName((PetscObject) *section, name);
228: #if 0
229:   SectionRealSetSection(*section, m->getRealSection(std::string(name)));
230:   SectionRealSetBundle(*section, m);
231: #endif
232:   return(0);
233: }

237: PetscErrorCode  DMSetFromOptions_Cartesian(DM dm)
238: {

243:   PetscOptionsHead("DMCartesian Options");
244:     /* Handle DMCartesian refinement */
245:     /* Handle associated vectors */
246:   PetscOptionsTail();
247:   return(0);
248: }

250: EXTERN_C_BEGIN
253: PetscErrorCode DMCreate_Cartesian(DM dm)
254: {
255:   DM_Cartesian  *mesh;

260:   PetscNewLog(dm, DM_Cartesian, &mesh);
261:   dm->data = mesh;

263:   new(&mesh->m) ALE::Obj<ALE::CartesianMesh>(PETSC_NULL);

265:   PetscStrallocpy(VECSTANDARD, &dm->vectype);
266:   dm->ops->globaltolocalbegin = 0;
267:   dm->ops->globaltolocalend   = 0;
268:   dm->ops->localtoglobalbegin = 0;
269:   dm->ops->localtoglobalend   = 0;
270:   dm->ops->createglobalvector = 0; /* DMCreateGlobalVector_Cartesian; */
271:   dm->ops->createlocalvector  = 0; /* DMCreateLocalVector_Cartesian; */
272:   dm->ops->createinterpolation   = DMCreateInterpolation_Cartesian;
273:   dm->ops->getcoloring        = 0;
274:   dm->ops->creatematrix          = 0; /* DMCreateMatrix_Cartesian; */
275:   dm->ops->refine             = DMRefine_Cartesian;
276:   dm->ops->coarsen            = DMCoarsen_Cartesian;
277:   dm->ops->refinehierarchy    = 0;
278:   dm->ops->coarsenhierarchy   = 0;
279:   dm->ops->getinjection       = 0;
280:   dm->ops->getaggregates      = 0;
281:   dm->ops->destroy            = DMDestroy_Cartesian;
282:   dm->ops->view               = DMView_Cartesian;
283:   dm->ops->setfromoptions     = DMSetFromOptions_Cartesian;
284:   dm->ops->setup              = 0;
285:   return(0);
286: }
287: EXTERN_C_END

291: /*@
292:   DMCartesianCreate - Creates a DMCartesian object.

294:   Collective on MPI_Comm

296:   Input Parameter:
297: . comm - The communicator for the DMCartesian object

299:   Output Parameter:
300: . mesh  - The DMCartesian object

302:   Level: beginner

304: .keywords: DMCartesian, create
305: @*/
306: PetscErrorCode DMCartesianCreate(MPI_Comm comm, DM *mesh)
307: {

312:   DMCreate(comm, mesh);
313:   DMSetType(*mesh, DMCARTESIAN);
314:   return(0);
315: }