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 createBoxMesh(MPI_Comm comm, const int dim, const bool structured, const bool interpolate, const int debug, DM *dm) {

 20:       if (structured) {
 21:         DA             da;
 22:         const PetscInt dof = 1;
 23:         const PetscInt pd  = PETSC_DECIDE;

 25:         if (dim == 2) {
 26:           DACreate2d(comm, DA_NONPERIODIC, DA_STENCIL_BOX, -3, -3, pd, pd, dof, 1, PETSC_NULL, PETSC_NULL, &da);
 27:         } else if (dim == 3) {
 28:           DACreate3d(comm, DA_NONPERIODIC, DA_STENCIL_BOX, -3, -3, -3, pd, pd, pd, dof, 1, PETSC_NULL, PETSC_NULL, PETSC_NULL, &da);
 29:         } else {
 30:           SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
 31:         }
 32:         DASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
 33:         *dm = (DM) da;
 34:       } else {
 35:         typedef PETSC_MESH_TYPE::point_type point_type;
 36:         ::Mesh mesh;
 37:         ::Mesh boundary;

 39:         MeshCreate(comm, &boundary);
 40:         Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
 41:         Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
 42:         std::map<point_type,point_type>  renumbering;
 43:         Obj<ALE::Mesh>                   mB;

 45:         meshBd->setSieve(sieve);
 46:         if (dim == 2) {
 47:           double lower[2] = {0.0, 0.0};
 48:           double upper[2] = {1.0, 1.0};
 49:           int    edges[2] = {2, 2};

 51:           mB = ALE::MeshBuilder<ALE::Mesh>::createSquareBoundary(comm, lower, upper, edges, debug);
 52:         } else if (dim == 3) {
 53:           double lower[3] = {0.0, 0.0, 0.0};
 54:           double upper[3] = {1.0, 1.0, 1.0};
 55:           int    faces[3] = {3, 3, 3};
 56: 
 57:           mB = ALE::MeshBuilder<ALE::Mesh>::createCubeBoundary(comm, lower, upper, faces, debug);
 58:         } else {
 59:           SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
 60:         }
 61:         ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
 62:         MeshSetMesh(boundary, meshBd);
 63:         MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
 64:         MeshDestroy(boundary);
 65:         *dm = (DM) mesh;
 66:       }
 67:       return(0);
 68:     };
 71:     static PetscErrorCode createReentrantBoxMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
 72:       typedef PETSC_MESH_TYPE::point_type point_type;
 73:       ::Mesh         mesh;
 74:       ::Mesh         boundary;

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

 84:       meshBd->setSieve(sieve);
 85:       if (dim == 2) {
 86:         double lower[2]  = {-1.0, -1.0};
 87:         double upper[2]  = {1.0, 1.0};
 88:         double offset[2] = {0.5, 0.5};

 90:         mB = ALE::MeshBuilder<ALE::Mesh>::createReentrantBoundary(comm, lower, upper, offset, debug);
 91:       } else if (dim == 3) {
 92:         double lower[3]  = {-1.0, -1.0, -1.0};
 93:         double upper[3]  = { 1.0,  1.0,  1.0};
 94:         double offset[3] = { 0.5,  0.5,  0.5};

 96:         mB = ALE::MeshBuilder<ALE::Mesh>::createFicheraCornerBoundary(comm, lower, upper, offset, debug);
 97:       } else {
 98:         SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
 99:       }
100:       ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
101:       MeshSetMesh(boundary, meshBd);
102:       MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
103:       MeshDestroy(boundary);
104:       *dm = (DM) mesh;
105:       return(0);
106:     };
109:     static PetscErrorCode createSphericalMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
110:       typedef PETSC_MESH_TYPE::point_type point_type;
111:       ::Mesh         mesh;
112:       ::Mesh         boundary;

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

122:       meshBd->setSieve(sieve);
123:       if (dim == 2) {
124:         mB = ALE::MeshBuilder<ALE::Mesh>::createCircularReentrantBoundary(comm, 100, 1.0, 1.0, debug);
125:       } else {
126:         SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
127:       }
128:       ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
129:       MeshSetMesh(boundary, meshBd);
130:       MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
131:       MeshDestroy(boundary);
132:       *dm = (DM) mesh;
133:       return(0);
134:     };
137:     static PetscErrorCode createReentrantSphericalMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
138:       typedef PETSC_MESH_TYPE::point_type point_type;
139:       ::Mesh         mesh;
140:       ::Mesh         boundary;

144:       MeshCreate(comm, &boundary);
145:       Obj<PETSC_MESH_TYPE>             meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
146:       Obj<PETSC_MESH_TYPE::sieve_type> sieve  = new PETSC_MESH_TYPE::sieve_type(comm, debug);
147:       std::map<point_type,point_type>  renumbering;
148:       Obj<ALE::Mesh>                   mB;

150:       meshBd->setSieve(sieve);
151:       if (dim == 2) {
152:         mB = ALE::MeshBuilder<ALE::Mesh>::createCircularReentrantBoundary(comm, 100, 1.0, 0.9, debug);
153:       } else {
154:         SETERRQ1(PETSC_ERR_SUP, "Dimension not supported: %d", dim);
155:       }
156:       ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
157:       MeshSetMesh(boundary, meshBd);
158:       MeshGenerate(boundary, (PetscTruth) interpolate, &mesh);
159:       MeshDestroy(boundary);
160:       *dm = (DM) mesh;
161:       return(0);
162:     };
165:     static PetscErrorCode MeshRefineSingularity(::Mesh mesh, double * singularity, double factor, ::Mesh *refinedMesh) {
166:       ALE::Obj<PETSC_MESH_TYPE> oldMesh;
167:       double              oldLimit;
168:       PetscErrorCode      ierr;

171:       MeshGetMesh(mesh, oldMesh);
172:       MeshCreate(oldMesh->comm(), refinedMesh);
173:       int dim = oldMesh->getDimension();
174:       oldLimit = oldMesh->getMaxVolume();
175:       //double oldLimInv = 1./oldLimit;
176:       double curLimit, tmpLimit;
177:       double minLimit = oldLimit/16384.;             //arbitrary;
178:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = oldMesh->getRealSection("coordinates");
179:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& volume_limits = oldMesh->getRealSection("volume_limits");
180:       volume_limits->setFiberDimension(oldMesh->heightStratum(0), 1);
181:       oldMesh->allocate(volume_limits);
182:       const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = oldMesh->heightStratum(0);
183:       PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin();
184:       PETSC_MESH_TYPE::label_sequence::iterator c_iter_end = cells->end();
185:       double centerCoords[dim];
186:       while (c_iter != c_iter_end) {
187:         const double * coords = oldMesh->restrictClosure(coordinates, *c_iter);
188:         for (int i = 0; i < dim; i++) {
189:           centerCoords[i] = 0;
190:           for (int j = 0; j < dim+1; j++) {
191:             centerCoords[i] += coords[j*dim+i];
192:           }
193:           centerCoords[i] = centerCoords[i]/(dim+1);
194:         }
195:         double dist = 0.;
196:         for (int i = 0; i < dim; i++) {
197:           dist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
198:         }
199:         if (dist > 0.) {
200:           dist = sqrt(dist);
201:           double mu = pow(dist, factor);
202:           //PetscPrintf(oldMesh->comm(), "%f\n", mu);
203:           tmpLimit = oldLimit*pow(mu, dim);
204:           if (tmpLimit > minLimit) {
205:             curLimit = tmpLimit;
206:           } else curLimit = minLimit;
207:         } else curLimit = minLimit;
208:         //PetscPrintf(oldMesh->comm(), "%f, %f\n", dist, tmpLimit);
209:         volume_limits->updatePoint(*c_iter, &curLimit);
210:         c_iter++;
211:       }
212: #ifdef PETSC_OPT_SIEVE
213:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, volume_limits, true);
214: #else
215:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMesh(oldMesh, volume_limits, true);
216: #endif
217:       MeshSetMesh(*refinedMesh, newMesh);
218:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& s = newMesh->getRealSection("default");
219:       const Obj<std::set<std::string> >& discs = oldMesh->getDiscretizations();

221:       for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
222:         newMesh->setDiscretization(*f_iter, oldMesh->getDiscretization(*f_iter));
223:       }
224:       newMesh->setupField(s);
225:       return(0);
226:     };
229:     static PetscErrorCode MeshRefineSingularity_Fichera(::Mesh mesh, double * singularity, double factor, ::Mesh *refinedMesh) {
230:       ALE::Obj<PETSC_MESH_TYPE> oldMesh;
231:       double              oldLimit;
232:       PetscErrorCode      ierr;

235:       MeshGetMesh(mesh, oldMesh);
236:       MeshCreate(oldMesh->comm(), refinedMesh);
237:       int dim = oldMesh->getDimension();
238:       oldLimit = oldMesh->getMaxVolume();
239:       //double oldLimInv = 1./oldLimit;
240:       double curLimit, tmpLimit;
241:       double minLimit = oldLimit/16384.;             //arbitrary;
242:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = oldMesh->getRealSection("coordinates");
243:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& volume_limits = oldMesh->getRealSection("volume_limits");
244:       volume_limits->setFiberDimension(oldMesh->heightStratum(0), 1);
245:       oldMesh->allocate(volume_limits);
246:       const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = oldMesh->heightStratum(0);
247:       PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin();
248:       PETSC_MESH_TYPE::label_sequence::iterator c_iter_end = cells->end();
249:       double centerCoords[dim];
250:       while (c_iter != c_iter_end) {
251:         const double * coords = oldMesh->restrictClosure(coordinates, *c_iter);
252:         for (int i = 0; i < dim; i++) {
253:           centerCoords[i] = 0;
254:           for (int j = 0; j < dim+1; j++) {
255:             centerCoords[i] += coords[j*dim+i];
256:           }
257:           centerCoords[i] = centerCoords[i]/(dim+1);
258:           //PetscPrintf(oldMesh->comm(), "%f, ", centerCoords[i]);
259:         }
260:         //PetscPrintf(oldMesh->comm(), "\n");
261:         double dist = 0.;
262:         double cornerdist = 0.;
263:         //HERE'S THE DIFFERENCE: if centercoords is less than the singularity coordinate for each direction, include that direction in the distance
264:         /*
265:           for (int i = 0; i < dim; i++) {
266:           if (centerCoords[i] <= singularity[i]) {
267:           dist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
268:           }
269:           }
270:         */
271:         //determine: the per-dimension distance: cases
272:         for (int i = 0; i < dim; i++) {
273:           cornerdist = 0.;
274:           if (centerCoords[i] > singularity[i]) {
275:             for (int j = 0; j < dim; j++) {
276:               if (j != i) cornerdist += (centerCoords[j] - singularity[j])*(centerCoords[j] - singularity[j]);
277:             }
278:             if (cornerdist < dist || dist == 0.) dist = cornerdist;
279:           }
280:         }
281:         //patch up AROUND the corner by minimizing between the distance from the relevant axis and the singular vertex
282:         double singdist = 0.;
283:         for (int i = 0; i < dim; i++) {
284:           singdist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
285:         }
286:         if (singdist < dist || dist == 0.) dist = singdist;
287:         if (dist > 0.) {
288:           dist = sqrt(dist);
289:           double mu = pow(dist, factor);
290:           //PetscPrintf(oldMesh->comm(), "%f, %f\n", mu, dist);
291:           tmpLimit = oldLimit*pow(mu, dim);
292:           if (tmpLimit > minLimit) {
293:             curLimit = tmpLimit;
294:           } else curLimit = minLimit;
295:         } else curLimit = minLimit;
296:         //PetscPrintf(oldMesh->comm(), "%f, %f\n", dist, tmpLimit);
297:         volume_limits->updatePoint(*c_iter, &curLimit);
298:         c_iter++;
299:       }
300: #ifdef PETSC_OPT_SIEVE
301:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, volume_limits, true);
302: #else
303:       ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMesh(oldMesh, volume_limits, true);
304: #endif
305:       MeshSetMesh(*refinedMesh, newMesh);
306:       const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& s = newMesh->getRealSection("default");
307:       const Obj<std::set<std::string> >& discs = oldMesh->getDiscretizations();

309:       for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
310:         newMesh->setDiscretization(*f_iter, oldMesh->getDiscretization(*f_iter));
311:       }
312:       newMesh->setupField(s);
313:       //  PetscPrintf(newMesh->comm(), "refined\n");
314:       return(0);
315:     };
316:   };
317: }

319: #endif