dm.c (97b5471de3aa6efc806ef4d184ef9c2a2ae62c92) dm.c (d0295fc027abbea29f13fa82c19cf92da8e9ba99)
1#include <petscvec.h>
1#include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
2#include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/
3#include <petsc/private/petscdsimpl.h> /*I "petscds.h" I*/
4#include <petscdmplex.h>
5#include <petscdmfield.h>
6#include <petscsf.h>
7#include <petscds.h>
8
9#if defined(PETSC_HAVE_VALGRIND)
10# include <valgrind/memcheck.h>
11#endif
12
13PetscClassId DM_CLASSID;
14PetscClassId DMLABEL_CLASSID;
15PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction, DM_CreateInjection, DM_CreateMatrix, DM_Load, DM_AdaptInterpolator;
16
17const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DMBoundaryType","DM_BOUNDARY_", NULL};
18const char *const DMBoundaryConditionTypes[] = {"INVALID","ESSENTIAL","NATURAL","INVALID","INVALID","ESSENTIAL_FIELD","NATURAL_FIELD","INVALID","INVALID","INVALID","NATURAL_RIEMANN","DMBoundaryConditionType","DM_BC_", NULL};
2#include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
3#include <petsc/private/dmlabelimpl.h> /*I "petscdmlabel.h" I*/
4#include <petsc/private/petscdsimpl.h> /*I "petscds.h" I*/
5#include <petscdmplex.h>
6#include <petscdmfield.h>
7#include <petscsf.h>
8#include <petscds.h>
9
10#if defined(PETSC_HAVE_VALGRIND)
11# include <valgrind/memcheck.h>
12#endif
13
14PetscClassId DM_CLASSID;
15PetscClassId DMLABEL_CLASSID;
16PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction, DM_CreateInjection, DM_CreateMatrix, DM_Load, DM_AdaptInterpolator;
17
18const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DMBoundaryType","DM_BOUNDARY_", NULL};
19const char *const DMBoundaryConditionTypes[] = {"INVALID","ESSENTIAL","NATURAL","INVALID","INVALID","ESSENTIAL_FIELD","NATURAL_FIELD","INVALID","INVALID","INVALID","NATURAL_RIEMANN","DMBoundaryConditionType","DM_BC_", NULL};
19const char *const DMPolytopeTypes[] = {"vertex", "segment", "tensor_segment", "triangle", "quadrilateral", "tensor_quad", "tetrahedron", "hexahedron", "triangular_prism", "tensor_triangular_prism", "tensor_quadrilateral_prism", "pyramid", "FV_ghost_cell", "interior_ghost_cell", "unknown", "invalid", "DMPolytopeType", "DM_POLYTOPE_", NULL};
20const char *const DMPolytopeTypes[] = {"vertex", "segment", "tensor_segment", "triangle", "quadrilateral", "tensor_quad", "tetrahedron", "hexahedron", "triangular_prism", "tensor_triangular_prism", "tensor_quadrilateral_prism", "FV_ghost_cell", "interior_ghost_cell", "unknown", "invalid", "DMPolytopeType", "DM_POLYTOPE_", NULL};
20/*@
21 DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
22
23 If you never call DMSetType() it will generate an
24 error when you try to use the vector.
25
26 Collective
27

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

1357 DM mdm;
1358
1359 ierr = MatGetDM(*mat,&mdm);CHKERRQ(ierr);
1360 if (!mdm) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"DM type '%s' did not attach the DM to the matrix\n",((PetscObject)dm)->type_name);
1361 }
1362 /* Handle nullspace and near nullspace */
1363 if (dm->Nf) {
1364 MatNullSpace nullSpace;
21/*@
22 DMCreate - Creates an empty DM object. The type can then be set with DMSetType().
23
24 If you never call DMSetType() it will generate an
25 error when you try to use the vector.
26
27 Collective
28

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

1358 DM mdm;
1359
1360 ierr = MatGetDM(*mat,&mdm);CHKERRQ(ierr);
1361 if (!mdm) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"DM type '%s' did not attach the DM to the matrix\n",((PetscObject)dm)->type_name);
1362 }
1363 /* Handle nullspace and near nullspace */
1364 if (dm->Nf) {
1365 MatNullSpace nullSpace;
1365 PetscInt Nf, f;
1366 PetscInt Nf;
1366
1367 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
1367
1368 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr);
1368 for (f = 0; f < Nf; ++f) {
1369 if (dm->nullspaceConstructors[f]) {
1370 ierr = (*dm->nullspaceConstructors[f])(dm, f, f, &nullSpace);CHKERRQ(ierr);
1369 if (Nf == 1) {
1370 if (dm->nullspaceConstructors[0]) {
1371 ierr = (*dm->nullspaceConstructors[0])(dm, 0, 0, &nullSpace);CHKERRQ(ierr);
1371 ierr = MatSetNullSpace(*mat, nullSpace);CHKERRQ(ierr);
1372 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
1372 ierr = MatSetNullSpace(*mat, nullSpace);CHKERRQ(ierr);
1373 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
1373 break;
1374 }
1374 }
1375 }
1376 for (f = 0; f < Nf; ++f) {
1377 if (dm->nearnullspaceConstructors[f]) {
1378 ierr = (*dm->nearnullspaceConstructors[f])(dm, f, f, &nullSpace);CHKERRQ(ierr);
1375 if (dm->nearnullspaceConstructors[0]) {
1376 ierr = (*dm->nearnullspaceConstructors[0])(dm, 0, 0, &nullSpace);CHKERRQ(ierr);
1379 ierr = MatSetNearNullSpace(*mat, nullSpace);CHKERRQ(ierr);
1380 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
1381 }
1382 }
1383 }
1384 ierr = PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);CHKERRQ(ierr);
1385 PetscFunctionReturn(0);
1386}

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

2487
2488@*/
2489PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2490{
2491 PetscSF sf;
2492 PetscErrorCode ierr;
2493 DMGlobalToLocalHookLink link;
2494
1377 ierr = MatSetNearNullSpace(*mat, nullSpace);CHKERRQ(ierr);
1378 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
1379 }
1380 }
1381 }
1382 ierr = PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);CHKERRQ(ierr);
1383 PetscFunctionReturn(0);
1384}

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

2485
2486@*/
2487PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2488{
2489 PetscSF sf;
2490 PetscErrorCode ierr;
2491 DMGlobalToLocalHookLink link;
2492
2493
2495 PetscFunctionBegin;
2496 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
2497 for (link=dm->gtolhook; link; link=link->next) {
2498 if (link->beginhook) {
2499 ierr = (*link->beginhook)(dm,g,mode,l,link->ctx);CHKERRQ(ierr);
2500 }
2501 }
2502 ierr = DMGetSectionSF(dm, &sf);CHKERRQ(ierr);
2503 if (sf) {
2504 const PetscScalar *gArray;
2505 PetscScalar *lArray;
2494 PetscFunctionBegin;
2495 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
2496 for (link=dm->gtolhook; link; link=link->next) {
2497 if (link->beginhook) {
2498 ierr = (*link->beginhook)(dm,g,mode,l,link->ctx);CHKERRQ(ierr);
2499 }
2500 }
2501 ierr = DMGetSectionSF(dm, &sf);CHKERRQ(ierr);
2502 if (sf) {
2503 const PetscScalar *gArray;
2504 PetscScalar *lArray;
2505 PetscMemType lmtype,gmtype;
2506
2507 if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2506
2507 if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2508 ierr = VecGetArrayInPlace(l, &lArray);CHKERRQ(ierr);
2509 ierr = VecGetArrayReadInPlace(g, &gArray);CHKERRQ(ierr);
2510 ierr = PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr);
2508 ierr = VecGetArrayInPlace_Internal(l, &lArray, &lmtype);CHKERRQ(ierr);
2509 ierr = VecGetArrayReadInPlace_Internal(g, &gArray, &gmtype);CHKERRQ(ierr);
2510 ierr = PetscSFBcastWithMemTypeBegin(sf, MPIU_SCALAR, gmtype, gArray, lmtype, lArray);CHKERRQ(ierr);
2511 ierr = VecRestoreArrayInPlace(l, &lArray);CHKERRQ(ierr);
2512 ierr = VecRestoreArrayReadInPlace(g, &gArray);CHKERRQ(ierr);
2513 } else {
2514 if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2515 ierr = (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr);
2516 }
2517 PetscFunctionReturn(0);
2518}

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

2536PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2537{
2538 PetscSF sf;
2539 PetscErrorCode ierr;
2540 const PetscScalar *gArray;
2541 PetscScalar *lArray;
2542 PetscBool transform;
2543 DMGlobalToLocalHookLink link;
2511 ierr = VecRestoreArrayInPlace(l, &lArray);CHKERRQ(ierr);
2512 ierr = VecRestoreArrayReadInPlace(g, &gArray);CHKERRQ(ierr);
2513 } else {
2514 if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2515 ierr = (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr);
2516 }
2517 PetscFunctionReturn(0);
2518}

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

2536PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2537{
2538 PetscSF sf;
2539 PetscErrorCode ierr;
2540 const PetscScalar *gArray;
2541 PetscScalar *lArray;
2542 PetscBool transform;
2543 DMGlobalToLocalHookLink link;
2544 PetscMemType lmtype,gmtype;
2544
2545 PetscFunctionBegin;
2546 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
2547 ierr = DMGetSectionSF(dm, &sf);CHKERRQ(ierr);
2548 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr);
2549 if (sf) {
2550 if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2551
2545
2546 PetscFunctionBegin;
2547 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
2548 ierr = DMGetSectionSF(dm, &sf);CHKERRQ(ierr);
2549 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr);
2550 if (sf) {
2551 if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2552
2552 ierr = VecGetArrayInPlace(l, &lArray);CHKERRQ(ierr);
2553 ierr = VecGetArrayReadInPlace(g, &gArray);CHKERRQ(ierr);
2553 ierr = VecGetArrayInPlace_Internal(l, &lArray, &lmtype);CHKERRQ(ierr);
2554 ierr = VecGetArrayReadInPlace_Internal(g, &gArray, &gmtype);CHKERRQ(ierr);
2554 ierr = PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr);
2555 ierr = VecRestoreArrayInPlace(l, &lArray);CHKERRQ(ierr);
2556 ierr = VecRestoreArrayReadInPlace(g, &gArray);CHKERRQ(ierr);
2557 if (transform) {ierr = DMPlexGlobalToLocalBasis(dm, l);CHKERRQ(ierr);}
2558 } else {
2559 if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2560 ierr = (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr);
2561 }

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

2713 PetscSF sf;
2714 PetscSection s, gs;
2715 DMLocalToGlobalHookLink link;
2716 Vec tmpl;
2717 const PetscScalar *lArray;
2718 PetscScalar *gArray;
2719 PetscBool isInsert, transform, l_inplace = PETSC_FALSE, g_inplace = PETSC_FALSE;
2720 PetscErrorCode ierr;
2555 ierr = PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr);
2556 ierr = VecRestoreArrayInPlace(l, &lArray);CHKERRQ(ierr);
2557 ierr = VecRestoreArrayReadInPlace(g, &gArray);CHKERRQ(ierr);
2558 if (transform) {ierr = DMPlexGlobalToLocalBasis(dm, l);CHKERRQ(ierr);}
2559 } else {
2560 if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2561 ierr = (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr);
2562 }

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

2714 PetscSF sf;
2715 PetscSection s, gs;
2716 DMLocalToGlobalHookLink link;
2717 Vec tmpl;
2718 const PetscScalar *lArray;
2719 PetscScalar *gArray;
2720 PetscBool isInsert, transform, l_inplace = PETSC_FALSE, g_inplace = PETSC_FALSE;
2721 PetscErrorCode ierr;
2722 PetscMemType lmtype=PETSC_MEMTYPE_HOST,gmtype=PETSC_MEMTYPE_HOST;
2721
2722 PetscFunctionBegin;
2723 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
2724 for (link=dm->ltoghook; link; link=link->next) {
2725 if (link->beginhook) {
2726 ierr = (*link->beginhook)(dm,l,mode,g,link->ctx);CHKERRQ(ierr);
2727 }
2728 }

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

2746 if (transform) {
2747 ierr = DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);CHKERRQ(ierr);
2748 ierr = VecCopy(l, tmpl);CHKERRQ(ierr);
2749 ierr = DMPlexLocalToGlobalBasis(dm, tmpl);CHKERRQ(ierr);
2750 ierr = VecGetArrayRead(tmpl, &lArray);CHKERRQ(ierr);
2751 } else if (isInsert) {
2752 ierr = VecGetArrayRead(l, &lArray);CHKERRQ(ierr);
2753 } else {
2723
2724 PetscFunctionBegin;
2725 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
2726 for (link=dm->ltoghook; link; link=link->next) {
2727 if (link->beginhook) {
2728 ierr = (*link->beginhook)(dm,l,mode,g,link->ctx);CHKERRQ(ierr);
2729 }
2730 }

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

2748 if (transform) {
2749 ierr = DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);CHKERRQ(ierr);
2750 ierr = VecCopy(l, tmpl);CHKERRQ(ierr);
2751 ierr = DMPlexLocalToGlobalBasis(dm, tmpl);CHKERRQ(ierr);
2752 ierr = VecGetArrayRead(tmpl, &lArray);CHKERRQ(ierr);
2753 } else if (isInsert) {
2754 ierr = VecGetArrayRead(l, &lArray);CHKERRQ(ierr);
2755 } else {
2754 ierr = VecGetArrayReadInPlace(l, &lArray);CHKERRQ(ierr);
2756 ierr = VecGetArrayReadInPlace_Internal(l, &lArray, &lmtype);CHKERRQ(ierr);
2755 l_inplace = PETSC_TRUE;
2756 }
2757 if (s && isInsert) {
2758 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr);
2759 } else {
2757 l_inplace = PETSC_TRUE;
2758 }
2759 if (s && isInsert) {
2760 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr);
2761 } else {
2760 ierr = VecGetArrayInPlace(g, &gArray);CHKERRQ(ierr);
2762 ierr = VecGetArrayInPlace_Internal(g, &gArray, &gmtype);CHKERRQ(ierr);
2761 g_inplace = PETSC_TRUE;
2762 }
2763 if (sf && !isInsert) {
2763 g_inplace = PETSC_TRUE;
2764 }
2765 if (sf && !isInsert) {
2764 ierr = PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);CHKERRQ(ierr);
2766 ierr = PetscSFReduceWithMemTypeBegin(sf, MPIU_SCALAR, lmtype, lArray, gmtype, gArray, MPIU_SUM);CHKERRQ(ierr);
2765 } else if (s && isInsert) {
2766 PetscInt gStart, pStart, pEnd, p;
2767
2768 ierr = DMGetGlobalSection(dm, &gs);CHKERRQ(ierr);
2769 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
2770 ierr = VecGetOwnershipRange(g, &gStart, NULL);CHKERRQ(ierr);
2771 for (p = pStart; p < pEnd; ++p) {
2772 PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

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

2859 PetscScalar *gArray;
2860 Vec tmpl;
2861
2862 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr);
2863 if (transform) {
2864 ierr = DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);CHKERRQ(ierr);
2865 ierr = VecGetArrayRead(tmpl, &lArray);CHKERRQ(ierr);
2866 } else {
2767 } else if (s && isInsert) {
2768 PetscInt gStart, pStart, pEnd, p;
2769
2770 ierr = DMGetGlobalSection(dm, &gs);CHKERRQ(ierr);
2771 ierr = PetscSectionGetChart(s, &pStart, &pEnd);CHKERRQ(ierr);
2772 ierr = VecGetOwnershipRange(g, &gStart, NULL);CHKERRQ(ierr);
2773 for (p = pStart; p < pEnd; ++p) {
2774 PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

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

2861 PetscScalar *gArray;
2862 Vec tmpl;
2863
2864 ierr = DMHasBasisTransform(dm, &transform);CHKERRQ(ierr);
2865 if (transform) {
2866 ierr = DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);CHKERRQ(ierr);
2867 ierr = VecGetArrayRead(tmpl, &lArray);CHKERRQ(ierr);
2868 } else {
2867 ierr = VecGetArrayReadInPlace(l, &lArray);CHKERRQ(ierr);
2869 ierr = VecGetArrayReadInPlace_Internal(l, &lArray, NULL);CHKERRQ(ierr);
2868 }
2870 }
2869 ierr = VecGetArrayInPlace(g, &gArray);CHKERRQ(ierr);
2871 ierr = VecGetArrayInPlace_Internal(g, &gArray, NULL);CHKERRQ(ierr);
2870 ierr = PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);CHKERRQ(ierr);
2871 if (transform) {
2872 ierr = VecRestoreArrayRead(tmpl, &lArray);CHKERRQ(ierr);
2873 ierr = DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);CHKERRQ(ierr);
2874 } else {
2875 ierr = VecRestoreArrayReadInPlace(l, &lArray);CHKERRQ(ierr);
2876 }
2877 ierr = VecRestoreArrayInPlace(g, &gArray);CHKERRQ(ierr);

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

6416 PetscFunctionBegin;
6417 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
6418 PetscValidHeaderSpecific(disc,PETSCFE_CLASSID,2);
6419
6420 ierr = DMGetCoordinateDM(dm, &cdmOld);CHKERRQ(ierr);
6421 /* Check current discretization is compatible */
6422 ierr = DMGetField(cdmOld, 0, NULL, &discOld);CHKERRQ(ierr);
6423 ierr = PetscObjectGetClassId(discOld, &classid);CHKERRQ(ierr);
2872 ierr = PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);CHKERRQ(ierr);
2873 if (transform) {
2874 ierr = VecRestoreArrayRead(tmpl, &lArray);CHKERRQ(ierr);
2875 ierr = DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);CHKERRQ(ierr);
2876 } else {
2877 ierr = VecRestoreArrayReadInPlace(l, &lArray);CHKERRQ(ierr);
2878 }
2879 ierr = VecRestoreArrayInPlace(g, &gArray);CHKERRQ(ierr);

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

6418 PetscFunctionBegin;
6419 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
6420 PetscValidHeaderSpecific(disc,PETSCFE_CLASSID,2);
6421
6422 ierr = DMGetCoordinateDM(dm, &cdmOld);CHKERRQ(ierr);
6423 /* Check current discretization is compatible */
6424 ierr = DMGetField(cdmOld, 0, NULL, &discOld);CHKERRQ(ierr);
6425 ierr = PetscObjectGetClassId(discOld, &classid);CHKERRQ(ierr);
6424 if (classid != PETSCFE_CLASSID) {
6425 if (classid == PETSC_CONTAINER_CLASSID) {
6426 PetscFE feLinear;
6427 DMPolytopeType ct;
6428 PetscInt dim, dE, cStart;
6429 PetscBool simplex;
6430
6431 /* Assume linear vertex coordinates */
6432 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
6433 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr);
6434 ierr = DMPlexGetHeightStratum(cdmOld, 0, &cStart, NULL);CHKERRQ(ierr);
6435 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr);
6436 switch (ct) {
6437 case DM_POLYTOPE_TRI_PRISM:
6438 case DM_POLYTOPE_TRI_PRISM_TENSOR:
6439 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot autoamtically create coordinate space for prisms");
6440 default: break;
6441 }
6442 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
6443 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, 1, -1, &feLinear);CHKERRQ(ierr);
6444 ierr = DMSetField(cdmOld, 0, NULL, (PetscObject) feLinear);CHKERRQ(ierr);
6445 ierr = PetscFEDestroy(&feLinear);CHKERRQ(ierr);
6446 ierr = DMCreateDS(cdmOld);CHKERRQ(ierr);
6447 } else {
6448 const char *discname;
6449
6450 ierr = PetscObjectGetType(discOld, &discname);CHKERRQ(ierr);
6451 SETERRQ1(PetscObjectComm(discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
6452 }
6453 }
6426 if (classid != PETSCFE_CLASSID) SETERRQ(PetscObjectComm(discOld), PETSC_ERR_SUP, "Discretization type not supported");
6454 /* Make a fresh clone of the coordinate DM */
6455 ierr = DMClone(cdmOld, &cdmNew);CHKERRQ(ierr);
6456 ierr = DMSetField(cdmNew, 0, NULL, (PetscObject) disc);CHKERRQ(ierr);
6457 ierr = DMCreateDS(cdmNew);CHKERRQ(ierr);
6458 /* Project the coordinate vector from old to new space */
6459 ierr = DMGetCoordinates(dm, &coordsOld);CHKERRQ(ierr);
6460 ierr = DMCreateGlobalVector(cdmNew, &coordsNew);CHKERRQ(ierr);
6461 ierr = DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);CHKERRQ(ierr);

--- 2671 unchanged lines hidden --- 		6427 /* Make a fresh clone of the coordinate DM */
6428 ierr = DMClone(cdmOld, &cdmNew);CHKERRQ(ierr);
6429 ierr = DMSetField(cdmNew, 0, NULL, (PetscObject) disc);CHKERRQ(ierr);
6430 ierr = DMCreateDS(cdmNew);CHKERRQ(ierr);
6431 /* Project the coordinate vector from old to new space */
6432 ierr = DMGetCoordinates(dm, &coordsOld);CHKERRQ(ierr);
6433 ierr = DMCreateGlobalVector(cdmNew, &coordsNew);CHKERRQ(ierr);
6434 ierr = DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);CHKERRQ(ierr);

--- 2671 unchanged lines hidden ---