Actual source code: DMBuilder.hh

  1: #ifndef included_ALE_DMBuilder_hh
  2: #define included_ALE_DMBuilder_hh

  4: #ifndef  included_ALE_Mesh_hh
  5: #include <Mesh.hh>
  6: #endif

  8: #include <petscmesh.hh>

 10: namespace ALE {

 12:   class DMBuilder {
 13:   public:
 16:     static PetscErrorCode createBasketMesh(MPI_Comm comm, const int dim, const bool structured, const bool interpolate, const int debug, DM *dm) {

 20:       if (structured) {
 21:         SETERRQ(PETSC_ERR_SUP, "Structured grids cannot handle boundary meshes");
 22:       } else {
 23:         typedef PETSC_MESH_TYPE::point_type point_type;
 24:         ::Mesh boundary;

 26:         MeshCreate(comm, &boundary);
 27:         Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
 28:         Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
 29:         std::map<point_type,point_type>  renumbering;
 30:         Obj<ALE::Mesh>                   mB;

 32:         meshBd->setSieve(sieve);
 33:         if (dim == 2) {
 34:           double lower[2] = {0.0, 0.0};
 35:           double upper[2] = {1.0, 1.0};
 36:           int    edges    = 2;

 38:           mB = ALE::MeshBuilder<ALE::Mesh>::createSquareBoundary(comm, lower, upper, edges, debug);
 39:         } else if (dim == 3) {
 40:           double lower[3] = {0.0, 0.0, 0.0};
 41:           double upper[3] = {1.0, 1.0, 1.0};
 42:           int    faces[3] = {3, 3, 3};
 43: 
 44:           mB = ALE::MeshBuilder<ALE::Mesh>::createCubeBoundary(comm, lower, upper, faces, debug);
 45:         } else {
 46:           SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
 47:         }
 48:         ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
 49:         MeshSetMesh(boundary, meshBd);
 50:         *dm = (DM) boundary;
 51:       }
 52:       return(0);
 53:     };
 56:     static PetscErrorCode createBoxMesh(MPI_Comm comm, const int dim, const bool structured, const bool interpolate, const int debug, DM *dm) {

 60:       if (structured) {
 61:         DA             da;
 62:         const PetscInt dof = 1;
 63:         const PetscInt pd  = PETSC_DECIDE;

 65:         if (dim == 2) {
 66:           DACreate2d(comm, DA_NONPERIODIC, DA_STENCIL_BOX, -3, -3, pd, pd, dof, 1, PETSC_NULL, PETSC_NULL, &da);
 67:         } else if (dim == 3) {
 68:           DACreate3d(comm, DA_NONPERIODIC, DA_STENCIL_BOX, -3, -3, -3, pd, pd, pd, dof, 1, PETSC_NULL, PETSC_NULL, PETSC_NULL, &da);
 69:         } else {
 70:           SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
 71:         }
 72:         DASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
 73:         *dm = (DM) da;
 74:       } else {
 75:         typedef PETSC_MESH_TYPE::point_type point_type;
 76:         ::Mesh mesh;
 77:         ::Mesh boundary;

 79:         MeshCreate(comm, &boundary);
 80:         Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
 81:         Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
 82:         std::map<point_type,point_type>  renumbering;
 83:         Obj<ALE::Mesh>                   mB;

 85:         meshBd->setSieve(sieve);
 86:         if (dim == 2) {
 87:           double lower[2] = {0.0, 0.0};
 88:           double upper[2] = {1.0, 1.0};
 89:           int    edges[2] = {2, 2};

 91:           mB = ALE::MeshBuilder<ALE::Mesh>::createSquareBoundary(comm, lower, upper, edges, debug);
 92:         } else if (dim == 3) {
 93:           double lower[3] = {0.0, 0.0, 0.0};
 94:           double upper[3] = {1.0, 1.0, 1.0};
 95:           int    faces[3] = {3, 3, 3};
 96: 
 97:           mB = ALE::MeshBuilder<ALE::Mesh>::createCubeBoundary(comm, lower, upper, faces, debug);
 98:         } else {
 99:           SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
100:         }
101:         ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
102:         MeshSetMesh(boundary, meshBd);
103:         MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
104:         MeshDestroy(boundary);
105:         *dm = (DM) mesh;
106:       }
107:       return(0);
108:     };
111:     static PetscErrorCode createReentrantBoxMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
112:       typedef PETSC_MESH_TYPE::point_type point_type;
113:       ::Mesh         mesh;
114:       ::Mesh         boundary;

118:       MeshCreate(comm, &boundary);
119:       Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
120:       Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
121:       std::map<point_type,point_type>  renumbering;
122:       Obj<ALE::Mesh>                   mB;

124:       meshBd->setSieve(sieve);
125:       if (dim == 2) {
126:         double lower[2]  = {-1.0, -1.0};
127:         double upper[2]  = {1.0, 1.0};
128:         double offset[2] = {0.5, 0.5};

130:         mB = ALE::MeshBuilder<ALE::Mesh>::createReentrantBoundary(comm, lower, upper, offset, debug);
131:       } else if (dim == 3) {
132:         double lower[3]  = {-1.0, -1.0, -1.0};
133:         double upper[3]  = { 1.0,  1.0,  1.0};
134:         double offset[3] = { 0.5,  0.5,  0.5};

136:         mB = ALE::MeshBuilder<ALE::Mesh>::createFicheraCornerBoundary(comm, lower, upper, offset, debug);
137:       } else {
138:         SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
139:       }
140:       ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
141:       MeshSetMesh(boundary, meshBd);
142:       MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
143:       MeshDestroy(boundary);
144:       *dm = (DM) mesh;
145:       return(0);
146:     };
149:     static PetscErrorCode createSphericalMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
150:       typedef PETSC_MESH_TYPE::point_type point_type;
151:       ::Mesh         mesh;
152:       ::Mesh         boundary;

156:       MeshCreate(comm, &boundary);
157:       Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
158:       Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
159:       std::map<point_type,point_type>  renumbering;
160:       Obj<ALE::Mesh>                   mB;

162:       meshBd->setSieve(sieve);
163:       if (dim == 2) {
164:         mB = ALE::MeshBuilder<ALE::Mesh>::createCircularReentrantBoundary(comm, 100, 1.0, 1.0, debug);
165:       } else {
166:         SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
167:       }
168:       ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
169:       MeshSetMesh(boundary, meshBd);
170:       MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
171:       MeshDestroy(boundary);
172:       *dm = (DM) mesh;
173:       return(0);
174:     };
177:     static PetscErrorCode createReentrantSphericalMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
178:       typedef PETSC_MESH_TYPE::point_type point_type;
179:       ::Mesh         mesh;
180:       ::Mesh         boundary;

184:       MeshCreate(comm, &boundary);
185:       Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
186:       Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
187:       std::map<point_type,point_type>  renumbering;
188:       Obj<ALE::Mesh>                   mB;

190:       meshBd->setSieve(sieve);
191:       if (dim == 2) {
192:         mB = ALE::MeshBuilder<ALE::Mesh>::createCircularReentrantBoundary(comm, 100, 1.0, 0.9, debug);
193:       } else {
194:         SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
195:       }
196:       ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
197:       MeshSetMesh(boundary, meshBd);
198:       MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
199:       MeshDestroy(boundary);
200:       *dm = (DM) mesh;
201:       return(0);
202:     };
205:     static PetscErrorCode MeshRefineSingularity(::Mesh mesh, double * singularity, double factor, ::Mesh *refinedMesh) {
206:       ALE::Obj<PETSC_MESH_TYPE> oldMesh;
207:       double              oldLimit;
208:       PetscErrorCode      ierr;

211:       MeshGetMesh(mesh, oldMesh);
212:       MeshCreate(oldMesh->comm(), refinedMesh);
213:       int dim = oldMesh->getDimension();
214:       oldLimit = oldMesh->getMaxVolume();
215:       //double oldLimInv = 1./oldLimit;
216:       double curLimit, tmpLimit;
217:       double minLimit = oldLimit/16384.;             //arbitrary;
218:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = oldMesh->getRealSection("coordinates");
219:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& volume_limits = oldMesh->getRealSection("volume_limits");
220:       volume_limits->setFiberDimension(oldMesh->heightStratum(0), 1);
221:       oldMesh->allocate(volume_limits);
222:       const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = oldMesh->heightStratum(0);
223:       PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin();
224:       PETSC_MESH_TYPE::label_sequence::iterator c_iter_end = cells->end();
225:       double centerCoords[dim];
226:       while (c_iter != c_iter_end) {
227:         const double * coords = oldMesh->restrictClosure(coordinates, *c_iter);
228:         for (int i = 0; i < dim; i++) {
229:           centerCoords[i] = 0;
230:           for (int j = 0; j < dim+1; j++) {
231:             centerCoords[i] += coords[j*dim+i];
232:           }
233:           centerCoords[i] = centerCoords[i]/(dim+1);
234:         }
235:         double dist = 0.;
236:         for (int i = 0; i < dim; i++) {
237:           dist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
238:         }
239:         if (dist > 0.) {
240:           dist = sqrt(dist);
241:           double mu = pow(dist, factor);
242:           //PetscPrintf(oldMesh->comm(), "%f\n", mu);
243:           tmpLimit = oldLimit*pow(mu, dim);
244:           if (tmpLimit > minLimit) {
245:             curLimit = tmpLimit;
246:           } else curLimit = minLimit;
247:         } else curLimit = minLimit;
248:         //PetscPrintf(oldMesh->comm(), "%f, %f\n", dist, tmpLimit);
249:         volume_limits->updatePoint(*c_iter, &curLimit);
250:         c_iter++;
251:       }
252: #ifdef PETSC_OPT_SIEVE
253:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, volume_limits, true);
254: #else
255:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMesh(oldMesh, volume_limits, true);
256: #endif
257:       MeshSetMesh(*refinedMesh, newMesh);
258:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& s = newMesh->getRealSection("default");
259:       const Obj<std::set<std::string> >& discs = oldMesh->getDiscretizations();

261:       for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
262:         newMesh->setDiscretization(*f_iter, oldMesh->getDiscretization(*f_iter));
263:       }
264:       newMesh->setupField(s);
265:       return(0);
266:     };
269:     static PetscErrorCode MeshRefineSingularity_Fichera(::Mesh mesh, double * singularity, double factor, ::Mesh *refinedMesh) {
270:       ALE::Obj<PETSC_MESH_TYPE> oldMesh;
271:       double              oldLimit;
272:       PetscErrorCode      ierr;

275:       MeshGetMesh(mesh, oldMesh);
276:       MeshCreate(oldMesh->comm(), refinedMesh);
277:       int dim = oldMesh->getDimension();
278:       oldLimit = oldMesh->getMaxVolume();
279:       //double oldLimInv = 1./oldLimit;
280:       double curLimit, tmpLimit;
281:       double minLimit = oldLimit/16384.;             //arbitrary;
282:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = oldMesh->getRealSection("coordinates");
283:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& volume_limits = oldMesh->getRealSection("volume_limits");
284:       volume_limits->setFiberDimension(oldMesh->heightStratum(0), 1);
285:       oldMesh->allocate(volume_limits);
286:       const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = oldMesh->heightStratum(0);
287:       PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin();
288:       PETSC_MESH_TYPE::label_sequence::iterator c_iter_end = cells->end();
289:       double centerCoords[dim];
290:       while (c_iter != c_iter_end) {
291:         const double * coords = oldMesh->restrictClosure(coordinates, *c_iter);
292:         for (int i = 0; i < dim; i++) {
293:           centerCoords[i] = 0;
294:           for (int j = 0; j < dim+1; j++) {
295:             centerCoords[i] += coords[j*dim+i];
296:           }
297:           centerCoords[i] = centerCoords[i]/(dim+1);
298:           //PetscPrintf(oldMesh->comm(), "%f, ", centerCoords[i]);
299:         }
300:         //PetscPrintf(oldMesh->comm(), "\n");
301:         double dist = 0.;
302:         double cornerdist = 0.;
303:         //HERE'S THE DIFFERENCE: if centercoords is less than the singularity coordinate for each direction, include that direction in the distance
304:         /*
305:           for (int i = 0; i < dim; i++) {
306:           if (centerCoords[i] <= singularity[i]) {
307:           dist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
308:           }
309:           }
310:         */
311:         //determine: the per-dimension distance: cases
312:         for (int i = 0; i < dim; i++) {
313:           cornerdist = 0.;
314:           if (centerCoords[i] > singularity[i]) {
315:             for (int j = 0; j < dim; j++) {
316:               if (j != i) cornerdist += (centerCoords[j] - singularity[j])*(centerCoords[j] - singularity[j]);
317:             }
318:             if (cornerdist < dist || dist == 0.) dist = cornerdist;
319:           }
320:         }
321:         //patch up AROUND the corner by minimizing between the distance from the relevant axis and the singular vertex
322:         double singdist = 0.;
323:         for (int i = 0; i < dim; i++) {
324:           singdist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
325:         }
326:         if (singdist < dist || dist == 0.) dist = singdist;
327:         if (dist > 0.) {
328:           dist = sqrt(dist);
329:           double mu = pow(dist, factor);
330:           //PetscPrintf(oldMesh->comm(), "%f, %f\n", mu, dist);
331:           tmpLimit = oldLimit*pow(mu, dim);
332:           if (tmpLimit > minLimit) {
333:             curLimit = tmpLimit;
334:           } else curLimit = minLimit;
335:         } else curLimit = minLimit;
336:         //PetscPrintf(oldMesh->comm(), "%f, %f\n", dist, tmpLimit);
337:         volume_limits->updatePoint(*c_iter, &curLimit);
338:         c_iter++;
339:       }
340: #ifdef PETSC_OPT_SIEVE
341:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, volume_limits, true);
342: #else
343:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMesh(oldMesh, volume_limits, true);
344: #endif
345:       MeshSetMesh(*refinedMesh, newMesh);
346:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& s = newMesh->getRealSection("default");
347:       const Obj<std::set<std::string> >& discs = oldMesh->getDiscretizations();

349:       for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
350:         newMesh->setDiscretization(*f_iter, oldMesh->getDiscretization(*f_iter));
351:       }
352:       newMesh->setupField(s);
353:       //  PetscPrintf(newMesh->comm(), "refined\n");
354:       return(0);
355:     };
356:   };
357: }

359: #endif