1 static char help[] = "Test metric utils in the uniform, isotropic case.\n\n"; 2 3 #include <petscdmplex.h> 4 5 static PetscErrorCode bowl(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 6 { 7 PetscInt d; 8 9 *u = 0.0; 10 for (d = 0; d < dim; d++) *u += 0.5*(x[d] - 0.5)*(x[d] - 0.5); 11 12 return 0; 13 } 14 15 int main(int argc, char **argv) { 16 DM dm, dmDist, dmAdapt; 17 DMLabel bdLabel = NULL; 18 DMPlexMetricCtx user; 19 MPI_Comm comm; 20 PetscBool uniform = PETSC_FALSE; 21 PetscErrorCode ierr; 22 PetscInt *faces, dim = 3, numEdges = 4, d; 23 PetscReal scaling = 1.0; 24 Vec metric; 25 26 /* Default parameter values */ 27 user.p = 1.0; 28 user.targetComplexity = 100.0; 29 user.isotropic = PETSC_FALSE; 30 31 /* Set up */ 32 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 33 comm = PETSC_COMM_WORLD; 34 ierr = PetscOptionsBegin(comm, "", "Mesh adaptation options", "DMPLEX");CHKERRQ(ierr); 35 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex60.c", dim, &dim, NULL, 2, 3);CHKERRQ(ierr); 36 ierr = PetscOptionsInt("-num_edges", "Number of edges on each boundary of the initial mesh", "ex60.c", numEdges, &numEdges, NULL);CHKERRQ(ierr); 37 ierr = PetscOptionsReal("-target", "Target metric complexity", "ex60.c", user.targetComplexity, &user.targetComplexity, NULL);CHKERRQ(ierr); 38 ierr = PetscOptionsReal("-norm_order", "Metric normalization order", "ex60.c", user.p, &user.p, NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsBool("-uniform", "Should the metric be assumed uniform?", "ex60.c", uniform, &uniform, NULL);CHKERRQ(ierr); 40 ierr = PetscOptionsBool("-isotropic", "Should the metric be assumed isotropic, or computed as a recovered Hessian?", "ex60.c", user.isotropic, &user.isotropic, NULL);CHKERRQ(ierr); 41 ierr = PetscOptionsEnd(); 42 43 /* Create box mesh */ 44 ierr = PetscMalloc1(dim, &faces);CHKERRQ(ierr); 45 for (d = 0; d < dim; ++d) faces[d] = numEdges; 46 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr); 47 ierr = PetscFree(faces);CHKERRQ(ierr); 48 49 /* Distribute mesh over processes */ 50 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr); 51 if (dmDist) { 52 ierr = DMDestroy(&dm);CHKERRQ(ierr); 53 dm = dmDist; 54 } 55 ierr = DMSetApplicationContext(dm, (void**)&user);CHKERRQ(ierr); 56 ierr = PetscObjectSetName((PetscObject) dm, "DM_init");CHKERRQ(ierr); 57 ierr = DMViewFromOptions(dm, NULL, "-initial_mesh_view");CHKERRQ(ierr); 58 59 /* Construct metric */ 60 if (uniform) { 61 if (user.isotropic) { ierr = DMPlexMetricCreateUniform(dm, 0, scaling, &metric);CHKERRQ(ierr); } 62 else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Uniform anisotropic metrics not supported."); 63 } 64 else { 65 DM dmIndi; 66 PetscFE fe; 67 Vec indicator; 68 69 /* Construct "error indicator" */ 70 ierr = DMClone(dm, &dmIndi);CHKERRQ(ierr); 71 ierr = PetscFECreateLagrange(comm, dim, 1, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 72 ierr = DMAddField(dmIndi, NULL, (PetscObject)fe);CHKERRQ(ierr); 73 ierr = DMCreateDS(dmIndi);CHKERRQ(ierr); 74 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 75 ierr = DMCreateLocalVector(dmIndi, &indicator);CHKERRQ(ierr); 76 if (user.isotropic) { 77 78 /* Isotropic case: just specify unity */ 79 ierr = VecSet(indicator, scaling);CHKERRQ(ierr); 80 ierr = DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric);CHKERRQ(ierr); 81 82 } else { 83 84 /* 'Anisotropic' case: approximate the identity by recovering the Hessian of a parabola */ 85 DM dmGrad; 86 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar*, void*) = {bowl}; 87 PetscInt pStart, pEnd, p; 88 PetscSection sec, secCoord; 89 Vec gradient; 90 91 /* Project the parabola into P1 space */ 92 ierr = DMProjectFunctionLocal(dmIndi, 0.0, funcs, NULL, INSERT_ALL_VALUES, indicator);CHKERRQ(ierr); 93 94 /* Approximate the gradient */ 95 ierr = DMClone(dmIndi, &dmGrad);CHKERRQ(ierr); 96 ierr = DMGetCoordinateSection(dmGrad, &secCoord);CHKERRQ(ierr); 97 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dmGrad), &sec);CHKERRQ(ierr); 98 ierr = PetscSectionSetNumFields(sec, 1);CHKERRQ(ierr); 99 ierr = PetscSectionSetFieldComponents(sec, 0, dim);CHKERRQ(ierr); 100 ierr = PetscSectionGetChart(secCoord, &pStart, &pEnd);CHKERRQ(ierr); 101 ierr = PetscSectionSetChart(sec, pStart, pEnd);CHKERRQ(ierr); 102 for (p = pStart; p < pEnd; ++p) { 103 ierr = PetscSectionSetDof(sec, p, dim);CHKERRQ(ierr); 104 ierr = PetscSectionSetFieldDof(sec, p, 0, dim);CHKERRQ(ierr); 105 } 106 ierr = PetscSectionSetUp(sec);CHKERRQ(ierr); 107 ierr = DMSetLocalSection(dmGrad, sec);CHKERRQ(ierr); 108 ierr = PetscSectionDestroy(&sec);CHKERRQ(ierr); 109 ierr = PetscFECreateLagrange(comm, dim, dim, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 110 ierr = DMClearFields(dmGrad);CHKERRQ(ierr); 111 ierr = DMAddField(dmGrad, NULL, (PetscObject)fe);CHKERRQ(ierr); 112 ierr = DMCreateDS(dmGrad);CHKERRQ(ierr); 113 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 114 ierr = DMCreateLocalVector(dmGrad, &gradient);CHKERRQ(ierr); 115 ierr = DMPlexComputeGradientClementInterpolant(dmIndi, indicator, gradient);CHKERRQ(ierr); 116 ierr = VecViewFromOptions(gradient, NULL, "-adapt_gradient_view");CHKERRQ(ierr); 117 118 /* Approximate the Hessian */ 119 ierr = DMGetCoordinateSection(dm, &secCoord);CHKERRQ(ierr); 120 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sec);CHKERRQ(ierr); 121 ierr = PetscSectionSetNumFields(sec, 1);CHKERRQ(ierr); 122 ierr = PetscSectionSetFieldComponents(sec, 0, dim*dim);CHKERRQ(ierr); 123 ierr = PetscSectionGetChart(secCoord, &pStart, &pEnd);CHKERRQ(ierr); 124 ierr = PetscSectionSetChart(sec, pStart, pEnd);CHKERRQ(ierr); 125 for (p = pStart; p < pEnd; ++p) { 126 ierr = PetscSectionSetDof(sec, p, dim*dim);CHKERRQ(ierr); 127 ierr = PetscSectionSetFieldDof(sec, p, 0, dim*dim);CHKERRQ(ierr); 128 } 129 ierr = PetscSectionSetUp(sec);CHKERRQ(ierr); 130 ierr = DMSetLocalSection(dm, sec);CHKERRQ(ierr); 131 ierr = PetscSectionDestroy(&sec);CHKERRQ(ierr); 132 ierr = PetscFECreateLagrange(comm, dim, dim*dim, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 133 ierr = DMClearFields(dm);CHKERRQ(ierr); 134 ierr = DMAddField(dm, NULL, (PetscObject)fe);CHKERRQ(ierr); 135 ierr = DMCreateDS(dm);CHKERRQ(ierr); 136 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 137 ierr = DMCreateLocalVector(dm, &metric);CHKERRQ(ierr); 138 ierr = DMPlexComputeGradientClementInterpolant(dmGrad, gradient, metric);CHKERRQ(ierr); 139 ierr = VecViewFromOptions(metric, NULL, "-adapt_hessian_view");CHKERRQ(ierr); 140 ierr = VecDestroy(&gradient);CHKERRQ(ierr); 141 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 142 } 143 ierr = VecDestroy(&indicator);CHKERRQ(ierr); 144 ierr = DMDestroy(&dmIndi);CHKERRQ(ierr); 145 } 146 147 /* Test metric routines */ 148 { 149 PetscReal errornorm, norm, tol = 1.0e-10, weights[2] = {0.8, 0.2}; 150 Vec metric1, metric2, metricComb; 151 Vec metrics[2]; 152 153 ierr = VecDuplicate(metric, &metric1);CHKERRQ(ierr); 154 ierr = VecSet(metric1, 0);CHKERRQ(ierr); 155 ierr = VecAXPY(metric1, 0.625, metric);CHKERRQ(ierr); 156 ierr = VecDuplicate(metric, &metric2);CHKERRQ(ierr); 157 ierr = VecSet(metric2, 0);CHKERRQ(ierr); 158 ierr = VecAXPY(metric2, 2.5, metric);CHKERRQ(ierr); 159 metrics[0] = metric1; 160 metrics[1] = metric2; 161 162 /* Test metric average */ 163 ierr = DMPlexMetricAverage(dm, 2, weights, metrics, &metricComb);CHKERRQ(ierr); 164 ierr = VecAXPY(metricComb, -1, metric);CHKERRQ(ierr); 165 ierr = VecNorm(metric, NORM_2, &norm);CHKERRQ(ierr); 166 ierr = VecNorm(metricComb, NORM_2, &errornorm);CHKERRQ(ierr); 167 errornorm /= norm; 168 if (errornorm > tol) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric average test failed (L2 error %f)", errornorm); 169 ierr = VecDestroy(&metricComb);CHKERRQ(ierr); 170 171 /* Test metric intersection */ 172 if (user.isotropic) { 173 ierr = DMPlexMetricIntersection(dm, 2, metrics, &metricComb);CHKERRQ(ierr); 174 ierr = VecAXPY(metricComb, -1, metric1);CHKERRQ(ierr); 175 ierr = VecNorm(metricComb, NORM_2, &errornorm);CHKERRQ(ierr); 176 errornorm /= norm; 177 if (errornorm > tol) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric intersection test failed (L2 error %f)", errornorm); 178 } 179 ierr = VecDestroy(&metric2);CHKERRQ(ierr); 180 ierr = VecDestroy(&metricComb);CHKERRQ(ierr); 181 ierr = VecCopy(metric, metric1);CHKERRQ(ierr); 182 183 /* Test metric SPD enforcement */ 184 ierr = DMPlexMetricEnforceSPD(dm, PETSC_TRUE, metric);CHKERRQ(ierr); 185 if (user.isotropic) { 186 ierr = VecAXPY(metric1, -1, metric);CHKERRQ(ierr); 187 ierr = VecNorm(metric1, NORM_2, &errornorm);CHKERRQ(ierr); 188 errornorm /= norm; 189 if (errornorm > tol) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric SPD enforcement test failed (L2 error %f)", errornorm); 190 } 191 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 192 193 /* Test metric normalization */ 194 ierr = DMPlexMetricNormalize(dm, metric, PETSC_TRUE, &metric1);CHKERRQ(ierr); 195 if (user.isotropic) { 196 scaling = PetscPowReal(user.targetComplexity, 2.0/dim); 197 ierr = DMPlexMetricCreateUniform(dm, 0, scaling, &metric2);CHKERRQ(ierr); 198 ierr = VecAXPY(metric2, -1, metric1);CHKERRQ(ierr); 199 ierr = VecNorm(metric2, NORM_2, &errornorm);CHKERRQ(ierr); 200 errornorm /= norm; 201 if (errornorm > tol) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric normalization test failed (L2 error %f)", errornorm); 202 } 203 ierr = VecCopy(metric1, metric);CHKERRQ(ierr); 204 ierr = VecDestroy(&metric2);CHKERRQ(ierr); 205 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 206 } 207 208 /* Adapt the mesh */ 209 ierr = DMAdaptMetric(dm, metric, bdLabel, &dmAdapt);CHKERRQ(ierr); 210 ierr = PetscObjectSetName((PetscObject) dmAdapt, "DM_adapted");CHKERRQ(ierr); 211 ierr = VecDestroy(&metric);CHKERRQ(ierr); 212 ierr = DMViewFromOptions(dmAdapt, NULL, "-adapted_mesh_view");CHKERRQ(ierr); 213 214 /* Compare DMs */ 215 ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 216 ierr = DMView(dmAdapt, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 217 218 /* Clean up */ 219 ierr = DMDestroy(&dmAdapt);CHKERRQ(ierr); 220 ierr = DMDestroy(&dm);CHKERRQ(ierr); 221 ierr = PetscFinalize(); 222 return 0; 223 } 224 225 /*TEST 226 227 build: 228 requires: pragmatic 229 230 test: 231 suffix: uniform_2d 232 args: -dim 2 -uniform -isotropic 233 test: 234 suffix: uniform_3d 235 args: -dim 3 -uniform -isotropic 236 test: 237 suffix: iso_2d 238 args: -dim 2 -isotropic 239 test: 240 suffix: iso_3d 241 args: -dim 3 -isotropic 242 test: 243 suffix: hessian_2d 244 args: -dim 2 245 test: 246 suffix: hessian_3d 247 args: -dim 3 248 249 TEST*/ 250