plexadapt.c (8713909f45e0191ddaf0544687829ca6f7972fa2) plexadapt.c (94ef8dde638caef1d0cd84a7dc8a2db65fcda8b6)
1#include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2#ifdef PETSC_HAVE_PRAGMATIC
3#include <pragmatic/cpragmatic.h>
4#endif
5
6/*
7 DMPlexRemesh_Internal - Generates a new mesh conforming to a metric field.
8
9 Input Parameters:
10+ dm - The DM object
1#include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2#ifdef PETSC_HAVE_PRAGMATIC
3#include <pragmatic/cpragmatic.h>
4#endif
5
6/*
7 DMPlexRemesh_Internal - Generates a new mesh conforming to a metric field.
8
9 Input Parameters:
10+ dm - The DM object
11. vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector
11. vertexMetric - The metric to which the mesh is adapted, defined vertex-wise.
12. remeshBd - Flag to allow boundary changes
13- bdLabelName - Label name for boundary tags which are preserved in dmNew, or NULL. Should not be "_boundary_".
14
15 Output Parameter:
16. dmNew - the new DM
17
18 Level: advanced
19
20.seealso: DMCoarsen(), DMRefine()
21*/
22PetscErrorCode DMPlexRemesh_Internal(DM dm, Vec vertexMetric, const char bdLabelName[], PetscBool remeshBd, DM *dmNew)
23{
12. remeshBd - Flag to allow boundary changes
13- bdLabelName - Label name for boundary tags which are preserved in dmNew, or NULL. Should not be "_boundary_".
14
15 Output Parameter:
16. dmNew - the new DM
17
18 Level: advanced
19
20.seealso: DMCoarsen(), DMRefine()
21*/
22PetscErrorCode DMPlexRemesh_Internal(DM dm, Vec vertexMetric, const char bdLabelName[], PetscBool remeshBd, DM *dmNew)
23{
24 MPI_Comm comm;
25 const char *bdName = "_boundary_";
24 const char *bdName = "_boundary_";
26 DM odm = dm, udm, cdm;
25 DM udm, cdm;
27 DMLabel bdLabel = NULL, bdLabelFull;
26 DMLabel bdLabel = NULL, bdLabelFull;
28 IS bdIS, globalVertexNum;
27 IS bdIS;
29 PetscSection coordSection;
30 Vec coordinates;
31 const PetscScalar *coords, *met;
28 PetscSection coordSection;
29 Vec coordinates;
30 const PetscScalar *coords, *met;
32 const PetscInt *bdFacesFull, *gV;
33 PetscInt *bdFaces, *bdFaceIds, *l2gv;
31 const PetscInt *bdFacesFull;
32 PetscInt *bdFaces, *bdFaceIds;
34 PetscReal *x, *y, *z, *metric;
35 PetscInt *cells;
33 PetscReal *x, *y, *z, *metric;
34 PetscInt *cells;
36 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
35 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, v;
37 PetscInt off, maxConeSize, numBdFaces, f, bdSize;
38 PetscBool flg;
39#ifdef PETSC_HAVE_PRAGMATIC
40 DMLabel bdLabelNew;
41 PetscScalar *coordsNew;
42 PetscInt *bdTags;
43 PetscReal *xNew[3] = {NULL, NULL, NULL};
44 PetscInt *cellsNew;
45 PetscInt d, numCellsNew, numVerticesNew;
46 PetscInt numCornersNew, fStart, fEnd;
47#endif
36 PetscInt off, maxConeSize, numBdFaces, f, bdSize;
37 PetscBool flg;
38#ifdef PETSC_HAVE_PRAGMATIC
39 DMLabel bdLabelNew;
40 PetscScalar *coordsNew;
41 PetscInt *bdTags;
42 PetscReal *xNew[3] = {NULL, NULL, NULL};
43 PetscInt *cellsNew;
44 PetscInt d, numCellsNew, numVerticesNew;
45 PetscInt numCornersNew, fStart, fEnd;
46#endif
48 PetscMPIInt numProcs;
49 PetscErrorCode ierr;
50
51 PetscFunctionBegin;
47 PetscErrorCode ierr;
48
49 PetscFunctionBegin;
52 /* Check for FEM adjacency flags */
53 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
54 PetscValidHeaderSpecific(vertexMetric, VEC_CLASSID, 2);
55 PetscValidCharPointer(bdLabelName, 3);
56 PetscValidPointer(dmNew, 5);
50 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
51 PetscValidHeaderSpecific(vertexMetric, VEC_CLASSID, 2);
52 PetscValidCharPointer(bdLabelName, 3);
53 PetscValidPointer(dmNew, 5);
57 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
58 ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr);
59 if (bdLabelName) {
60 size_t len;
61
62 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr);
54 if (bdLabelName) {
55 size_t len;
56
57 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr);
63 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
58 if (flg) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
64 ierr = PetscStrlen(bdLabelName, &len);CHKERRQ(ierr);
65 if (len) {
66 ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr);
59 ierr = PetscStrlen(bdLabelName, &len);CHKERRQ(ierr);
60 if (len) {
61 ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr);
67 if (!bdLabel) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Label \"%s\" does not exist in DM", bdLabelName);
62 if (!bdLabel) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Label \"%s\" does not exist in DM", bdLabelName);
68 }
69 }
63 }
64 }
70 /* Add overlap for Pragmatic */
71#if 0
72 /* Check for overlap by looking for cell in the SF */
73 if (!overlapped) {
74 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr);
75 if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);}
76 }
77#endif
78 /* Get mesh information */
79 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
80 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
81 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
82 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
83 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr);
84 numCells = cEnd - cStart;
85 numVertices = vEnd - vStart;
86 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr);
87 for (c = 0, coff = 0; c < numCells; ++c) {
88 const PetscInt *cone;
89 PetscInt coneSize, cl;
90
91 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
92 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
93 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
94 }
65 /* Get mesh information */
66 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
67 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
68 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
69 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr);
70 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr);
71 numCells = cEnd - cStart;
72 numVertices = vEnd - vStart;
73 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr);
74 for (c = 0, coff = 0; c < numCells; ++c) {
75 const PetscInt *cone;
76 PetscInt coneSize, cl;
77
78 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr);
79 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr);
80 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart;
81 }
95 ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr);
96 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr);
97 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr);
98 for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
99 if (gV[v] >= 0) ++numLocVertices;
100 l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v];
101 }
102 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr);
103 ierr = DMDestroy(&udm);CHKERRQ(ierr);
104 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
105 ierr = DMGetDefaultSection(cdm, &coordSection);CHKERRQ(ierr);
106 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
107 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr);
108 for (v = vStart; v < vEnd; ++v) {
109 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
110 x[v-vStart] = PetscRealPart(coords[off+0]);

--- 31 unchanged lines hidden (view full) ---

142 else {bdFaceIds[f] = 1;}
143 }
144 ierr = ISDestroy(&bdIS);CHKERRQ(ierr);
145 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr);
146 /* Get metric */
147 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr);
148 for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
149 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr);
82 ierr = DMDestroy(&udm);CHKERRQ(ierr);
83 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
84 ierr = DMGetDefaultSection(cdm, &coordSection);CHKERRQ(ierr);
85 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
86 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr);
87 for (v = vStart; v < vEnd; ++v) {
88 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
89 x[v-vStart] = PetscRealPart(coords[off+0]);

--- 31 unchanged lines hidden (view full) ---

121 else {bdFaceIds[f] = 1;}
122 }
123 ierr = ISDestroy(&bdIS);CHKERRQ(ierr);
124 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr);
125 /* Get metric */
126 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr);
127 for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]);
128 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr);
150#if 0
151 /* Destroy overlap mesh */
152 ierr = DMDestroy(&dm);CHKERRQ(ierr);
153#endif
154 /* Create new mesh */
155#ifdef PETSC_HAVE_PRAGMATIC
156 switch (dim) {
157 case 2:
129 /* Create new mesh */
130#ifdef PETSC_HAVE_PRAGMATIC
131 switch (dim) {
132 case 2:
158 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break;
133 pragmatic_2d_init(&numVertices, &numCells, cells, x, y);break;
159 case 3:
134 case 3:
160 pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break;
161 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
135 pragmatic_3d_init(&numVertices, &numCells, cells, x, y, z);break;
136 default: SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
162 }
163 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
164 pragmatic_set_metric(metric);
165 pragmatic_adapt(remeshBd ? 1 : 0);
137 }
138 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds);
139 pragmatic_set_metric(metric);
140 pragmatic_adapt(remeshBd ? 1 : 0);
166 ierr = PetscFree(l2gv);CHKERRQ(ierr);
167 /* Read out mesh */
141 /* Read out mesh */
168 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew);
142 pragmatic_get_info(&numVerticesNew, &numCellsNew);
169 ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr);
170 switch (dim) {
171 case 2:
172 numCornersNew = 3;
173 ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr);
143 ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr);
144 switch (dim) {
145 case 2:
146 numCornersNew = 3;
147 ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr);
174 pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]);
148 pragmatic_get_coords_2d(xNew[0], xNew[1]);
175 break;
176 case 3:
177 numCornersNew = 4;
178 ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr);
149 break;
150 case 3:
151 numCornersNew = 4;
152 ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr);
179 pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]);
153 pragmatic_get_coords_3d(xNew[0], xNew[1], xNew[2]);
180 break;
181 default:
154 break;
155 default:
182 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
156 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim);
183 }
184 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];}
185 ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr);
186 pragmatic_get_elements(cellsNew);
157 }
158 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];}
159 ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr);
160 pragmatic_get_elements(cellsNew);
187 ierr = DMPlexCreateFromCellListParallel(comm, dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);CHKERRQ(ierr);
161 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, dmNew);CHKERRQ(ierr);
188 /* Read out boundary label */
189 pragmatic_get_boundaryTags(&bdTags);
190 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr);
191 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr);
192 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr);
193 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr);
194 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr);
195 for (c = cStart; c < cEnd; ++c) {
196 /* Only for simplicial meshes */
197 coff = (c-cStart)*(dim+1);
162 /* Read out boundary label */
163 pragmatic_get_boundaryTags(&bdTags);
164 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr);
165 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr);
166 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr);
167 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr);
168 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr);
169 for (c = cStart; c < cEnd; ++c) {
170 /* Only for simplicial meshes */
171 coff = (c-cStart)*(dim+1);
198 /* d is the local cell number of the vertex opposite to the face we are marking */
199 for (d = 0; d < dim+1; ++d) {
200 if (bdTags[coff+d]) {
172 for (d = 0; d < dim+1; ++d) {
173 if (bdTags[coff+d]) {
201 const PetscInt perm[4][4] = {{-1, -1, -1, -1}, {-1, -1, -1, -1}, {1, 2, 0, -1}, {3, 2, 1, 0}}; /* perm[d] = face opposite */
202 const PetscInt *cone;
174 const PetscInt *faces;
175 PetscInt numFaces, vertices[3], p;
203
176
204 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
205 ierr = DMPlexGetCone(*dmNew, c, &cone);CHKERRQ(ierr);
206 ierr = DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);CHKERRQ(ierr);
177 /* Mark face opposite to this vertex */
178 for (p = 0, v = 0; p < dim+1; ++p) if (p != d) {vertices[v] = cellsNew[coff+p] + vStart; ++v;}
179 ierr = DMPlexGetFullJoin(*dmNew, dim, vertices, &numFaces, &faces);CHKERRQ(ierr);
180 if (numFaces != 1) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot find boundary face in new cell %D for vertex %D(%D) with tag %D", c, cellsNew[coff+d], d, bdTags[coff+d]);
181 if (faces[0] < fStart || faces[0] >= fEnd) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Boundary point %D in new cell %D for vertex %D(%D) with tag %D is not a face", faces[0], c, cellsNew[coff+d], d, bdTags[coff+d]);
182 ierr = DMLabelSetValue(bdLabelNew, faces[0], bdTags[coff+d]);CHKERRQ(ierr);
183 ierr = DMPlexRestoreJoin(*dmNew, dim, vertices, &numFaces, &faces);CHKERRQ(ierr);
207 }
208 }
209 }
210 /* Cleanup */
211 switch (dim) {
212 case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);break;
213 case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);break;
214 }
215 ierr = PetscFree(cellsNew);CHKERRQ(ierr);
216 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr);
217 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
218 ierr = PetscFree(coordsNew);CHKERRQ(ierr);
219 pragmatic_finalize();
220#else
184 }
185 }
186 }
187 /* Cleanup */
188 switch (dim) {
189 case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);break;
190 case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);break;
191 }
192 ierr = PetscFree(cellsNew);CHKERRQ(ierr);
193 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr);
194 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr);
195 ierr = PetscFree(coordsNew);CHKERRQ(ierr);
196 pragmatic_finalize();
197#else
221 SETERRQ(comm, PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
198 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic.");
222#endif
223 PetscFunctionReturn(0);
224}
225
199#endif
200 PetscFunctionReturn(0);
201}
202
226/*@
203/*@C
227 DMPlexAdapt - Generates a mesh adapted to the specified metric field using the pragmatic library.
228
229 Input Parameters:
230+ dm - The DM object
231. metric - The metric to which the mesh is adapted, defined vertex-wise.
232- bdyLabelName - Label name for boundary tags. These will be preserved in the output mesh. bdyLabelName should be "" (empty string) if there is no such label, and should be different from "boundary".
233
234 Output Parameter:

--- 19 unchanged lines hidden ---
204 DMPlexAdapt - Generates a mesh adapted to the specified metric field using the pragmatic library.
205
206 Input Parameters:
207+ dm - The DM object
208. metric - The metric to which the mesh is adapted, defined vertex-wise.
209- bdyLabelName - Label name for boundary tags. These will be preserved in the output mesh. bdyLabelName should be "" (empty string) if there is no such label, and should be different from "boundary".
210
211 Output Parameter:

--- 19 unchanged lines hidden ---