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