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, rgLabel = NULL; 18 MPI_Comm comm; 19 PetscBool uniform = PETSC_FALSE, isotropic = PETSC_FALSE; 20 PetscErrorCode ierr; 21 PetscInt *faces, dim = 3, d; 22 PetscReal scaling = 1.0; 23 Vec metric; 24 25 /* Set up */ 26 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 27 comm = PETSC_COMM_WORLD; 28 ierr = PetscOptionsBegin(comm, "", "Mesh adaptation options", "DMPLEX");CHKERRQ(ierr); 29 ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex60.c", dim, &dim, NULL, 2, 3);CHKERRQ(ierr); 30 ierr = PetscOptionsBool("-uniform", "Should the metric be assumed uniform?", "ex60.c", uniform, &uniform, NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsBool("-isotropic", "Should the metric be assumed isotropic, or computed as a recovered Hessian?", "ex60.c", isotropic, &isotropic, NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsEnd(); 33 34 /* Create box mesh */ 35 ierr = PetscMalloc1(dim, &faces);CHKERRQ(ierr); 36 for (d = 0; d < dim; ++d) faces[d] = 4; 37 ierr = DMPlexCreateBoxMesh(comm, dim, PETSC_TRUE, faces, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr); 38 ierr = PetscFree(faces);CHKERRQ(ierr); 39 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 40 41 /* Distribute mesh over processes */ 42 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr); 43 if (dmDist) { 44 ierr = DMDestroy(&dm);CHKERRQ(ierr); 45 dm = dmDist; 46 } 47 ierr = PetscObjectSetName((PetscObject) dm, "DM_init");CHKERRQ(ierr); 48 ierr = DMViewFromOptions(dm, NULL, "-initial_mesh_view");CHKERRQ(ierr); 49 50 /* Construct metric */ 51 if (uniform) { 52 if (isotropic) { ierr = DMPlexMetricCreateUniform(dm, 0, scaling, &metric);CHKERRQ(ierr); } 53 else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Uniform anisotropic metrics not supported."); 54 } 55 else { 56 DM dmIndi; 57 PetscFE fe; 58 Vec indicator; 59 60 /* Construct "error indicator" */ 61 ierr = DMClone(dm, &dmIndi);CHKERRQ(ierr); 62 ierr = PetscFECreateLagrange(comm, dim, 1, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 63 ierr = DMSetField(dmIndi, 0, NULL, (PetscObject)fe);CHKERRQ(ierr); 64 ierr = DMCreateDS(dmIndi);CHKERRQ(ierr); 65 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 66 ierr = DMCreateLocalVector(dmIndi, &indicator);CHKERRQ(ierr); 67 if (isotropic) { 68 69 /* Isotropic case: just specify unity */ 70 ierr = VecSet(indicator, scaling);CHKERRQ(ierr); 71 ierr = DMPlexMetricCreateIsotropic(dm, 0, indicator, &metric);CHKERRQ(ierr); 72 73 } else { 74 75 /* 'Anisotropic' case: approximate the identity by recovering the Hessian of a parabola */ 76 DM dmGrad; 77 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar*, void*) = {bowl}; 78 Vec gradient; 79 80 /* Project the parabola into P1 space */ 81 ierr = DMProjectFunctionLocal(dmIndi, 0.0, funcs, NULL, INSERT_ALL_VALUES, indicator);CHKERRQ(ierr); 82 83 /* Approximate the gradient */ 84 ierr = DMClone(dmIndi, &dmGrad);CHKERRQ(ierr); 85 ierr = PetscFECreateLagrange(comm, dim, dim, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);CHKERRQ(ierr); 86 ierr = DMSetField(dmGrad, 0, NULL, (PetscObject)fe);CHKERRQ(ierr); 87 ierr = DMCreateDS(dmGrad);CHKERRQ(ierr); 88 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 89 ierr = DMCreateLocalVector(dmGrad, &gradient);CHKERRQ(ierr); 90 ierr = DMPlexComputeGradientClementInterpolant(dmIndi, indicator, gradient);CHKERRQ(ierr); 91 ierr = VecViewFromOptions(gradient, NULL, "-adapt_gradient_view");CHKERRQ(ierr); 92 93 /* Approximate the Hessian */ 94 ierr = DMPlexMetricCreate(dm, 0, &metric);CHKERRQ(ierr); 95 ierr = DMPlexComputeGradientClementInterpolant(dmGrad, gradient, metric);CHKERRQ(ierr); 96 ierr = VecViewFromOptions(metric, NULL, "-adapt_hessian_view");CHKERRQ(ierr); 97 ierr = VecDestroy(&gradient);CHKERRQ(ierr); 98 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 99 } 100 ierr = VecDestroy(&indicator);CHKERRQ(ierr); 101 ierr = DMDestroy(&dmIndi);CHKERRQ(ierr); 102 } 103 104 /* Test metric routines */ 105 { 106 DM dmDet; 107 PetscReal errornorm, norm, tol = 1.0e-10, weights[2] = {0.8, 0.2}; 108 Vec metric1, metric2, metricComb, determinant; 109 Vec metrics[2]; 110 111 ierr = VecDuplicate(metric, &metric1);CHKERRQ(ierr); 112 ierr = VecSet(metric1, 0);CHKERRQ(ierr); 113 ierr = VecAXPY(metric1, 0.625, metric);CHKERRQ(ierr); 114 ierr = VecDuplicate(metric, &metric2);CHKERRQ(ierr); 115 ierr = VecSet(metric2, 0);CHKERRQ(ierr); 116 ierr = VecAXPY(metric2, 2.5, metric);CHKERRQ(ierr); 117 metrics[0] = metric1; 118 metrics[1] = metric2; 119 120 /* Test metric average */ 121 ierr = DMPlexMetricAverage(dm, 2, weights, metrics, &metricComb);CHKERRQ(ierr); 122 ierr = VecAXPY(metricComb, -1, metric);CHKERRQ(ierr); 123 ierr = VecNorm(metric, NORM_2, &norm);CHKERRQ(ierr); 124 ierr = VecNorm(metricComb, NORM_2, &errornorm);CHKERRQ(ierr); 125 errornorm /= norm; 126 ierr = PetscPrintf(comm, "Metric average L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 127 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric average test failed"); 128 ierr = VecDestroy(&metricComb);CHKERRQ(ierr); 129 130 /* Test metric intersection */ 131 if (isotropic) { 132 ierr = DMPlexMetricIntersection(dm, 2, metrics, &metricComb);CHKERRQ(ierr); 133 ierr = VecAXPY(metricComb, -1, metric1);CHKERRQ(ierr); 134 ierr = VecNorm(metricComb, NORM_2, &errornorm);CHKERRQ(ierr); 135 errornorm /= norm; 136 ierr = PetscPrintf(comm, "Metric intersection L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 137 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric intersection test failed"); 138 } 139 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 140 ierr = VecDestroy(&metric2);CHKERRQ(ierr); 141 ierr = VecDestroy(&metricComb);CHKERRQ(ierr); 142 143 /* Test metric SPD enforcement */ 144 ierr = DMPlexMetricEnforceSPD(dm, metric, PETSC_TRUE, PETSC_TRUE, &metric1, &determinant);CHKERRQ(ierr); 145 if (isotropic) { 146 Vec err; 147 148 ierr = VecDuplicate(determinant, &err);CHKERRQ(ierr); 149 ierr = VecSet(err, 1.0);CHKERRQ(ierr); 150 ierr = VecNorm(err, NORM_2, &norm);CHKERRQ(ierr); 151 ierr = VecAXPY(err, -1, determinant);CHKERRQ(ierr); 152 ierr = VecNorm(err, NORM_2, &errornorm);CHKERRQ(ierr); 153 ierr = VecDestroy(&err);CHKERRQ(ierr); 154 errornorm /= norm; 155 ierr = PetscPrintf(comm, "Metric determinant L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 156 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Determinant is not unit"); 157 ierr = VecAXPY(metric1, -1, metric);CHKERRQ(ierr); 158 ierr = VecNorm(metric1, NORM_2, &errornorm);CHKERRQ(ierr); 159 errornorm /= norm; 160 ierr = PetscPrintf(comm, "Metric SPD enforcement L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 161 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric SPD enforcement test failed"); 162 } 163 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 164 ierr = VecGetDM(determinant, &dmDet);CHKERRQ(ierr); 165 ierr = VecDestroy(&determinant);CHKERRQ(ierr); 166 ierr = DMDestroy(&dmDet);CHKERRQ(ierr); 167 168 /* Test metric normalization */ 169 ierr = DMPlexMetricNormalize(dm, metric, PETSC_TRUE, PETSC_TRUE, &metric1);CHKERRQ(ierr); 170 if (isotropic) { 171 PetscReal target; 172 173 ierr = DMPlexMetricGetTargetComplexity(dm, &target);CHKERRQ(ierr); 174 scaling = PetscPowReal(target, 2.0/dim); 175 ierr = DMPlexMetricCreateUniform(dm, 0, scaling, &metric2);CHKERRQ(ierr); 176 ierr = VecAXPY(metric2, -1, metric1);CHKERRQ(ierr); 177 ierr = VecNorm(metric2, NORM_2, &errornorm);CHKERRQ(ierr); 178 errornorm /= norm; 179 ierr = PetscPrintf(comm, "Metric normalization L2 error: %.4f%%\n", 100*errornorm);CHKERRQ(ierr); 180 if (errornorm > tol) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Metric normalization test failed"); 181 } 182 ierr = VecCopy(metric1, metric);CHKERRQ(ierr); 183 ierr = VecDestroy(&metric2);CHKERRQ(ierr); 184 ierr = VecDestroy(&metric1);CHKERRQ(ierr); 185 } 186 187 /* Adapt the mesh */ 188 ierr = DMAdaptMetric(dm, metric, bdLabel, rgLabel, &dmAdapt);CHKERRQ(ierr); 189 ierr = DMDestroy(&dm);CHKERRQ(ierr); 190 ierr = PetscObjectSetName((PetscObject) dmAdapt, "DM_adapted");CHKERRQ(ierr); 191 ierr = VecDestroy(&metric);CHKERRQ(ierr); 192 ierr = DMViewFromOptions(dmAdapt, NULL, "-adapted_mesh_view");CHKERRQ(ierr); 193 194 /* Clean up */ 195 ierr = DMDestroy(&dmAdapt);CHKERRQ(ierr); 196 ierr = PetscFinalize(); 197 return 0; 198 } 199 200 /*TEST 201 202 test: 203 suffix: uniform_2d_pragmatic 204 requires: pragmatic 205 args: -dm_plex_metric_target_complexity 100 -dim 2 -dm_adaptor pragmatic -uniform -isotropic 206 test: 207 suffix: uniform_3d_pragmatic 208 requires: pragmatic tetgen 209 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor pragmatic -uniform -isotropic 210 test: 211 suffix: iso_2d_pragmatic 212 requires: pragmatic 213 args: -dm_plex_metric_target_complexity 100 -dim 2 -dm_adaptor pragmatic -isotropic 214 test: 215 suffix: iso_3d_pragmatic 216 requires: pragmatic tetgen 217 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor pragmatic -isotropic 218 test: 219 suffix: hessian_2d_pragmatic 220 requires: pragmatic 221 args: -dm_plex_metric_target_complexity 100 -dim 2 -dm_adaptor pragmatic 222 test: 223 suffix: hessian_3d_pragmatic 224 requires: pragmatic tetgen 225 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor pragmatic 226 test: 227 suffix: uniform_2d_mmg 228 requires: mmg 229 args: -dm_plex_metric_target_complexity 100 -dim 2 -dm_adaptor mmg -uniform -isotropic 230 test: 231 suffix: uniform_3d_mmg 232 requires: mmg tetgen 233 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor mmg -uniform -isotropic 234 test: 235 suffix: iso_2d_mmg 236 requires: mmg 237 args: -dm_plex_metric_target_complexity 100 -dim 2 -dm_adaptor mmg -isotropic 238 test: 239 suffix: iso_3d_mmg 240 requires: mmg tetgen 241 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor mmg -isotropic 242 test: 243 suffix: hessian_2d_mmg 244 requires: mmg 245 args: -dm_plex_metric_target_complexity 100 -dim 2 -dm_adaptor mmg 246 test: 247 suffix: hessian_3d_mmg 248 requires: mmg tetgen 249 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor mmg 250 test: 251 suffix: uniform_3d_parmmg 252 requires: parmmg tetgen 253 nsize: 2 254 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor parmmg -uniform -isotropic 255 test: 256 suffix: iso_3d_parmmg 257 requires: parmmg tetgen 258 nsize: 2 259 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor parmmg -isotropic 260 test: 261 suffix: hessian_3d_parmmg 262 requires: parmmg tetgen 263 nsize: 2 264 args: -dm_plex_metric_target_complexity 100 -dim 3 -dm_adaptor parmmg 265 266 TEST*/ 267