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