Actual source code: DMBuilder.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_DMBuilder_hh
2: #define included_ALE_DMBuilder_hh
4: #ifndef included_ALE_Mesh_hh
5: #include <sieve/Mesh.hh>
6: #endif
8: #include <petscdmmesh.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) {
17: typedef PETSC_MESH_TYPE::real_section_type::value_type real;
21: if (structured) {
22: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Structured grids cannot handle boundary meshes");
23: } else {
24: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
25: typedef PETSC_MESH_TYPE::point_type point_type;
26: DM boundary;
28: DMCreate(comm, &boundary);
29: DMSetType(boundary, DMMESH);
30: Obj<PETSC_MESH_TYPE> meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
31: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, debug);
32: std::map<point_type,point_type> renumbering;
33: Obj<FlexMesh> mB;
35: meshBd->setSieve(sieve);
36: if (dim == 2) {
37: real lower[2] = {0.0, 0.0};
38: real upper[2] = {1.0, 1.0};
39: int edges = 2;
41: mB = ALE::MeshBuilder<FlexMesh>::createSquareBoundary(comm, lower, upper, edges, debug);
42: } else if (dim == 3) {
43: real lower[3] = {0.0, 0.0, 0.0};
44: real upper[3] = {1.0, 1.0, 1.0};
45: int faces[3] = {3, 3, 3};
47: mB = ALE::MeshBuilder<FlexMesh>::createCubeBoundary(comm, lower, upper, faces, debug);
48: } else {
49: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "Dimension not supported: %d", dim);
50: }
51: ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
52: DMMeshSetMesh(boundary, meshBd);
53: *dm = boundary;
54: }
55: return(0);
56: };
59: static PetscErrorCode createBoxMesh(MPI_Comm comm, const int dim, const bool structured, const bool interpolate, const int debug, DM *dm) {
60: typedef PETSC_MESH_TYPE::real_section_type::value_type real;
64: if (structured) {
65: DM da;
66: const PetscInt dof = 1;
67: const PetscInt pd = PETSC_DECIDE;
69: if (dim == 2) {
70: DMDACreate2d(comm, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, -3, -3, pd, pd, dof, 1, PETSC_NULL, PETSC_NULL, &da);
71: } else if (dim == 3) {
72: DMDACreate3d(comm, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX, -3, -3, -3, pd, pd, pd, dof, 1, PETSC_NULL, PETSC_NULL, PETSC_NULL, &da);
73: } else {
74: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "Dimension not supported: %d", dim);
75: }
76: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
77: *dm = da;
78: } else {
79: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
80: typedef PETSC_MESH_TYPE::point_type point_type;
81: DM mesh;
82: DM boundary;
84: DMCreate(comm, &boundary);
85: DMSetType(boundary, DMMESH);
86: Obj<PETSC_MESH_TYPE> meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
87: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, debug);
88: std::map<point_type,point_type> renumbering;
89: Obj<FlexMesh> mB;
91: meshBd->setSieve(sieve);
92: if (dim == 2) {
93: real lower[2] = {0.0, 0.0};
94: real upper[2] = {1.0, 1.0};
95: int edges[2] = {2, 2};
97: mB = ALE::MeshBuilder<FlexMesh>::createSquareBoundary(comm, lower, upper, edges, debug);
98: } else if (dim == 3) {
99: real lower[3] = {0.0, 0.0, 0.0};
100: real upper[3] = {1.0, 1.0, 1.0};
101: int faces[3] = {3, 3, 3};
103: mB = ALE::MeshBuilder<FlexMesh>::createCubeBoundary(comm, lower, upper, faces, debug);
104: } else {
105: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "Dimension not supported: %d", dim);
106: }
107: ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
108: DMMeshSetMesh(boundary, meshBd);
109: DMMeshGenerate(boundary, (PetscBool) interpolate, &mesh);
110: DMDestroy(&boundary);
111: *dm = mesh;
112: }
113: return(0);
114: };
117: static PetscErrorCode createReentrantBoxMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
118: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
119: typedef PETSC_MESH_TYPE::point_type point_type;
120: typedef PETSC_MESH_TYPE::real_section_type::value_type real;
121: DM mesh;
122: DM boundary;
126: DMCreate(comm, &boundary);
127: DMSetType(boundary, DMMESH);
128: Obj<PETSC_MESH_TYPE> meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
129: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, debug);
130: std::map<point_type,point_type> renumbering;
131: Obj<FlexMesh> mB;
133: meshBd->setSieve(sieve);
134: if (dim == 2) {
135: real lower[2] = {-1.0, -1.0};
136: real upper[2] = {1.0, 1.0};
137: real offset[2] = {0.5, 0.5};
139: mB = ALE::MeshBuilder<FlexMesh>::createReentrantBoundary(comm, lower, upper, offset, debug);
140: } else if (dim == 3) {
141: real lower[3] = {-1.0, -1.0, -1.0};
142: real upper[3] = { 1.0, 1.0, 1.0};
143: real offset[3] = { 0.5, 0.5, 0.5};
145: mB = ALE::MeshBuilder<FlexMesh>::createFicheraCornerBoundary(comm, lower, upper, offset, debug);
146: } else {
147: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "Dimension not supported: %d", dim);
148: }
149: ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
150: DMMeshSetMesh(boundary, meshBd);
151: DMMeshGenerate(boundary, (PetscBool) interpolate, &mesh);
152: DMDestroy(&boundary);
153: *dm = mesh;
154: return(0);
155: };
158: static PetscErrorCode createSphericalMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
159: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
160: typedef PETSC_MESH_TYPE::point_type point_type;
161: DM mesh;
162: DM boundary;
166: DMCreate(comm, &boundary);
167: DMSetType(boundary, DMMESH);
168: Obj<PETSC_MESH_TYPE> meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
169: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, debug);
170: std::map<point_type,point_type> renumbering;
171: Obj<FlexMesh> mB;
173: meshBd->setSieve(sieve);
174: if (dim == 2) {
175: mB = ALE::MeshBuilder<FlexMesh>::createCircularReentrantBoundary(comm, 100, 1.0, 1.0, debug);
176: } else {
177: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "Dimension not supported: %d", dim);
178: }
179: ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
180: DMMeshSetMesh(boundary, meshBd);
181: DMMeshGenerate(boundary, (PetscBool) interpolate, &mesh);
182: DMDestroy(&boundary);
183: *dm = mesh;
184: return(0);
185: };
188: static PetscErrorCode createReentrantSphericalMesh(MPI_Comm comm, const int dim, const bool interpolate, const int debug, DM *dm) {
189: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
190: typedef PETSC_MESH_TYPE::point_type point_type;
191: DM mesh;
192: DM boundary;
196: DMCreate(comm, &boundary);
197: DMSetType(boundary, DMMESH);
198: Obj<PETSC_MESH_TYPE> meshBd = new PETSC_MESH_TYPE(comm, dim-1, debug);
199: Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(comm, debug);
200: std::map<point_type,point_type> renumbering;
201: Obj<FlexMesh> mB;
203: meshBd->setSieve(sieve);
204: if (dim == 2) {
205: mB = ALE::MeshBuilder<FlexMesh>::createCircularReentrantBoundary(comm, 100, 1.0, 0.9, debug);
206: } else {
207: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP, "Dimension not supported: %d", dim);
208: }
209: ALE::ISieveConverter::convertMesh(*mB, *meshBd, renumbering, false);
210: DMMeshSetMesh(boundary, meshBd);
211: DMMeshGenerate(boundary, (PetscBool) interpolate, &mesh);
212: DMDestroy(&boundary);
213: *dm = mesh;
214: return(0);
215: };
218: static PetscErrorCode MeshRefineSingularity(DM mesh, double * singularity, double factor, DM *refinedMesh) {
219: typedef PETSC_MESH_TYPE::real_section_type::value_type real;
220: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
221: double oldLimit;
222: PetscErrorCode ierr;
225: DMMeshGetMesh(mesh, oldMesh);
226: DMCreate(oldMesh->comm(), refinedMesh);
227: DMSetType(*refinedMesh, DMMESH);
228: int dim = oldMesh->getDimension();
229: oldLimit = oldMesh->getMaxVolume();
230: //double oldLimInv = 1./oldLimit;
231: real curLimit, tmpLimit;
232: real minLimit = oldLimit/16384.; //arbitrary;
233: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = oldMesh->getRealSection("coordinates");
234: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& volume_limits = oldMesh->getRealSection("volume_limits");
235: volume_limits->setFiberDimension(oldMesh->heightStratum(0), 1);
236: oldMesh->allocate(volume_limits);
237: const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = oldMesh->heightStratum(0);
238: PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin();
239: PETSC_MESH_TYPE::label_sequence::iterator c_iter_end = cells->end();
240: real centerCoords[dim];
241: while (c_iter != c_iter_end) {
242: const real * coords = oldMesh->restrictClosure(coordinates, *c_iter);
243: for (int i = 0; i < dim; i++) {
244: centerCoords[i] = 0;
245: for (int j = 0; j < dim+1; j++) {
246: centerCoords[i] += coords[j*dim+i];
247: }
248: centerCoords[i] = centerCoords[i]/(dim+1);
249: }
250: real dist = 0.;
251: for (int i = 0; i < dim; i++) {
252: dist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
253: }
254: if (dist > 0.) {
255: dist = sqrt(dist);
256: real mu = pow(dist, factor);
257: //PetscPrintf(oldMesh->comm(), "%f\n", mu);
258: tmpLimit = oldLimit*pow(mu, dim);
259: if (tmpLimit > minLimit) {
260: curLimit = tmpLimit;
261: } else curLimit = minLimit;
262: } else curLimit = minLimit;
263: //PetscPrintf(oldMesh->comm(), "%f, %f\n", dist, tmpLimit);
264: volume_limits->updatePoint(*c_iter, &curLimit);
265: c_iter++;
266: }
267: #ifdef PETSC_OPT_SIEVE
268: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, volume_limits, true);
269: #else
270: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMesh(oldMesh, volume_limits, true);
271: #endif
272: DMMeshSetMesh(*refinedMesh, newMesh);
273: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& s = newMesh->getRealSection("default");
274: const Obj<std::set<std::string> >& discs = oldMesh->getDiscretizations();
276: for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
277: newMesh->setDiscretization(*f_iter, oldMesh->getDiscretization(*f_iter));
278: }
279: newMesh->setupField(s);
280: return(0);
281: };
284: static PetscErrorCode MeshRefineSingularity_Fichera(DM mesh, double * singularity, double factor, DM *refinedMesh) {
285: typedef PETSC_MESH_TYPE::real_section_type::value_type real;
286: ALE::Obj<PETSC_MESH_TYPE> oldMesh;
287: real oldLimit;
288: PetscErrorCode ierr;
291: DMMeshGetMesh(mesh, oldMesh);
292: DMCreate(oldMesh->comm(), refinedMesh);
293: DMSetType(*refinedMesh, DMMESH);
294: int dim = oldMesh->getDimension();
295: oldLimit = oldMesh->getMaxVolume();
296: //double oldLimInv = 1./oldLimit;
297: real curLimit, tmpLimit;
298: real minLimit = oldLimit/16384.; //arbitrary;
299: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& coordinates = oldMesh->getRealSection("coordinates");
300: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& volume_limits = oldMesh->getRealSection("volume_limits");
301: volume_limits->setFiberDimension(oldMesh->heightStratum(0), 1);
302: oldMesh->allocate(volume_limits);
303: const ALE::Obj<PETSC_MESH_TYPE::label_sequence>& cells = oldMesh->heightStratum(0);
304: PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin();
305: PETSC_MESH_TYPE::label_sequence::iterator c_iter_end = cells->end();
306: real centerCoords[dim];
307: while (c_iter != c_iter_end) {
308: const real *coords = oldMesh->restrictClosure(coordinates, *c_iter);
309: for (int i = 0; i < dim; i++) {
310: centerCoords[i] = 0;
311: for (int j = 0; j < dim+1; j++) {
312: centerCoords[i] += coords[j*dim+i];
313: }
314: centerCoords[i] = centerCoords[i]/(dim+1);
315: //PetscPrintf(oldMesh->comm(), "%f, ", centerCoords[i]);
316: }
317: //PetscPrintf(oldMesh->comm(), "\n");
318: real dist = 0.;
319: real cornerdist = 0.;
320: //HERE'S THE DIFFERENCE: if centercoords is less than the singularity coordinate for each direction, include that direction in the distance
321: /*
322: for (int i = 0; i < dim; i++) {
323: if (centerCoords[i] <= singularity[i]) {
324: dist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
325: }
326: }
327: */
328: //determine: the per-dimension distance: cases
329: for (int i = 0; i < dim; i++) {
330: cornerdist = 0.;
331: if (centerCoords[i] > singularity[i]) {
332: for (int j = 0; j < dim; j++) {
333: if (j != i) cornerdist += (centerCoords[j] - singularity[j])*(centerCoords[j] - singularity[j]);
334: }
335: if (cornerdist < dist || dist == 0.) dist = cornerdist;
336: }
337: }
338: //patch up AROUND the corner by minimizing between the distance from the relevant axis and the singular vertex
339: real singdist = 0.;
340: for (int i = 0; i < dim; i++) {
341: singdist += (centerCoords[i] - singularity[i])*(centerCoords[i] - singularity[i]);
342: }
343: if (singdist < dist || dist == 0.) dist = singdist;
344: if (dist > 0.) {
345: dist = sqrt(dist);
346: real mu = pow(dist, factor);
347: //PetscPrintf(oldMesh->comm(), "%f, %f\n", mu, dist);
348: tmpLimit = oldLimit*pow(mu, dim);
349: if (tmpLimit > minLimit) {
350: curLimit = tmpLimit;
351: } else curLimit = minLimit;
352: } else curLimit = minLimit;
353: //PetscPrintf(oldMesh->comm(), "%f, %f\n", dist, tmpLimit);
354: volume_limits->updatePoint(*c_iter, &curLimit);
355: c_iter++;
356: }
357: #ifdef PETSC_OPT_SIEVE
358: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMeshV(oldMesh, volume_limits, true);
359: #else
360: ALE::Obj<PETSC_MESH_TYPE> newMesh = ALE::Generator<PETSC_MESH_TYPE>::refineMesh(oldMesh, volume_limits, true);
361: #endif
362: DMMeshSetMesh(*refinedMesh, newMesh);
363: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& s = newMesh->getRealSection("default");
364: const Obj<std::set<std::string> >& discs = oldMesh->getDiscretizations();
366: for(std::set<std::string>::const_iterator f_iter = discs->begin(); f_iter != discs->end(); ++f_iter) {
367: newMesh->setDiscretization(*f_iter, oldMesh->getDiscretization(*f_iter));
368: }
369: newMesh->setupField(s);
370: // PetscPrintf(newMesh->comm(), "refined\n");
371: return(0);
372: };
373: };
374: }
376: #endif