plex.c (0629556046e1c17b069549e1fb8598edbd177d43) plex.c (d02c7345ebed081476ce6bcdfd784b34c024d6ef)
1#include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2#include <petsc/private/dmlabelimpl.h>
3#include <petsc/private/isimpl.h>
4#include <petsc/private/vecimpl.h>
5#include <petsc/private/glvisvecimpl.h>
6#include <petscsf.h>
7#include <petscds.h>
8#include <petscdraw.h>

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

2677 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMInterpolateSolution_C", NULL));
2678 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertTimeDerivativeBoundaryValues_C", NULL));
2679 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", NULL));
2680 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexDistributeGetDefault_C", NULL));
2681 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexDistributeSetDefault_C", NULL));
2682 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", NULL));
2683 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderGetDefault_C", NULL));
2684 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderSetDefault_C", NULL));
1#include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2#include <petsc/private/dmlabelimpl.h>
3#include <petsc/private/isimpl.h>
4#include <petsc/private/vecimpl.h>
5#include <petsc/private/glvisvecimpl.h>
6#include <petscsf.h>
7#include <petscds.h>
8#include <petscdraw.h>

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

2677 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMInterpolateSolution_C", NULL));
2678 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexInsertTimeDerivativeBoundaryValues_C", NULL));
2679 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", NULL));
2680 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexDistributeGetDefault_C", NULL));
2681 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexDistributeSetDefault_C", NULL));
2682 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", NULL));
2683 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderGetDefault_C", NULL));
2684 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderSetDefault_C", NULL));
2685 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderSectionGetDefault_C", NULL));
2686 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexReorderSectionSetDefault_C", NULL));
2685 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", NULL));
2686 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexSetOverlap_C", NULL));
2687 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetUseCeed_C", NULL));
2688 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexSetUseCeed_C", NULL));
2689 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", NULL));
2690 if (--mesh->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2691 PetscCall(PetscSectionDestroy(&mesh->coneSection));
2692 PetscCall(PetscFree(mesh->cones));

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

2721 if (mesh->metricCtx) PetscCall(PetscFree(mesh->metricCtx));
2722 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
2723 PetscCall(PetscFree(mesh));
2724 PetscFunctionReturn(PETSC_SUCCESS);
2725}
2726
2727PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
2728{
2687 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetOverlap_C", NULL));
2688 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexSetOverlap_C", NULL));
2689 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexGetUseCeed_C", NULL));
2690 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMPlexSetUseCeed_C", NULL));
2691 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "DMGetIsoperiodicPointSF_C", NULL));
2692 if (--mesh->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2693 PetscCall(PetscSectionDestroy(&mesh->coneSection));
2694 PetscCall(PetscFree(mesh->cones));

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

2723 if (mesh->metricCtx) PetscCall(PetscFree(mesh->metricCtx));
2724 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
2725 PetscCall(PetscFree(mesh));
2726 PetscFunctionReturn(PETSC_SUCCESS);
2727}
2728
2729PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J)
2730{
2729 PetscSection sectionGlobal;
2731 PetscSection sectionGlobal, sectionLocal;
2730 PetscInt bs = -1, mbs;
2731 PetscInt localSize, localStart = 0;
2732 PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isMatIS;
2733 MatType mtype;
2734 ISLocalToGlobalMapping ltog;
2735
2736 PetscFunctionBegin;
2737 PetscCall(MatInitializePackage());
2738 mtype = dm->mattype;
2732 PetscInt bs = -1, mbs;
2733 PetscInt localSize, localStart = 0;
2734 PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isMatIS;
2735 MatType mtype;
2736 ISLocalToGlobalMapping ltog;
2737
2738 PetscFunctionBegin;
2739 PetscCall(MatInitializePackage());
2740 mtype = dm->mattype;
2741 PetscCall(DMGetLocalSection(dm, &sectionLocal));
2739 PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
2740 /* PetscCall(PetscSectionGetStorageSize(sectionGlobal, &localSize)); */
2741 PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize));
2742 PetscCallMPI(MPI_Exscan(&localSize, &localStart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
2743 PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
2744 PetscCall(MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE));
2745 PetscCall(MatSetType(*J, mtype));
2746 PetscCall(MatSetFromOptions(*J));

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

2767 for (p = pStart; p < pEnd; ++p) {
2768 switch (dm->blocking_type) {
2769 case DM_BLOCKING_TOPOLOGICAL_POINT: { // One block per topological point
2770 PetscInt bdof, offset;
2771
2772 PetscCall(PetscSectionGetDof(sectionGlobal, p, &dof));
2773 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &offset));
2774 PetscCall(PetscSectionGetConstraintDof(sectionGlobal, p, &cdof));
2742 PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
2743 /* PetscCall(PetscSectionGetStorageSize(sectionGlobal, &localSize)); */
2744 PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize));
2745 PetscCallMPI(MPI_Exscan(&localSize, &localStart, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
2746 PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
2747 PetscCall(MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE));
2748 PetscCall(MatSetType(*J, mtype));
2749 PetscCall(MatSetFromOptions(*J));

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

2770 for (p = pStart; p < pEnd; ++p) {
2771 switch (dm->blocking_type) {
2772 case DM_BLOCKING_TOPOLOGICAL_POINT: { // One block per topological point
2773 PetscInt bdof, offset;
2774
2775 PetscCall(PetscSectionGetDof(sectionGlobal, p, &dof));
2776 PetscCall(PetscSectionGetOffset(sectionGlobal, p, &offset));
2777 PetscCall(PetscSectionGetConstraintDof(sectionGlobal, p, &cdof));
2775 for (PetscInt i = 0; i < dof - cdof; i++) pblocks[offset - localStart + i] = dof - cdof;
2778 for (PetscInt i = 0; i < dof - cdof; ++i) pblocks[offset - localStart + i] = dof - cdof;
2779 // Signal block concatenation
2780 if (dof - cdof && sectionLocal->blockStarts && !PetscBTLookup(sectionLocal->blockStarts, p)) pblocks[offset - localStart] = -(dof - cdof);
2776 dof = dof < 0 ? -(dof + 1) : dof;
2777 bdof = cdof && (dof - cdof) ? 1 : dof;
2778 if (dof) {
2779 if (bs < 0) {
2780 bs = bdof;
2781 } else if (bs != bdof) {
2782 bs = 1;
2783 }

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

2822 PetscCall(PetscCalloc4(localSize / bs, &dnz, localSize / bs, &onz, localSize / bs, &dnzu, localSize / bs, &onzu));
2823 PetscCall(DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix));
2824 PetscCall(PetscFree4(dnz, onz, dnzu, onzu));
2825 }
2826 { // Consolidate blocks
2827 PetscInt nblocks = 0;
2828 for (PetscInt i = 0; i < localSize; i += PetscMax(1, pblocks[i])) {
2829 if (pblocks[i] == 0) continue;
2781 dof = dof < 0 ? -(dof + 1) : dof;
2782 bdof = cdof && (dof - cdof) ? 1 : dof;
2783 if (dof) {
2784 if (bs < 0) {
2785 bs = bdof;
2786 } else if (bs != bdof) {
2787 bs = 1;
2788 }

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

2827 PetscCall(PetscCalloc4(localSize / bs, &dnz, localSize / bs, &onz, localSize / bs, &dnzu, localSize / bs, &onzu));
2828 PetscCall(DMPlexPreallocateOperator(dm, bs, dnz, onz, dnzu, onzu, *J, fillMatrix));
2829 PetscCall(PetscFree4(dnz, onz, dnzu, onzu));
2830 }
2831 { // Consolidate blocks
2832 PetscInt nblocks = 0;
2833 for (PetscInt i = 0; i < localSize; i += PetscMax(1, pblocks[i])) {
2834 if (pblocks[i] == 0) continue;
2830 pblocks[nblocks++] = pblocks[i]; // nblocks always <= i
2835 // Negative block size indicates the blocks should be concatenated
2836 if (pblocks[i] < 0) {
2837 pblocks[i] = -pblocks[i];
2838 pblocks[nblocks - 1] += pblocks[i];
2839 } else {
2840 pblocks[nblocks++] = pblocks[i]; // nblocks always <= i
2841 }
2831 for (PetscInt j = 1; j < pblocks[i]; j++) PetscCheck(pblocks[i + j] == pblocks[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Block of size %" PetscInt_FMT " mismatches entry %" PetscInt_FMT, pblocks[i], pblocks[i + j]);
2832 }
2833 PetscCall(MatSetVariableBlockSizes(*J, nblocks, pblocks));
2834 }
2835 PetscCall(PetscFree(pblocks));
2836 }
2837 PetscCall(MatSetDM(*J, dm));
2838 PetscFunctionReturn(PETSC_SUCCESS);

--- 7695 unchanged lines hidden --- 		2842 for (PetscInt j = 1; j < pblocks[i]; j++) PetscCheck(pblocks[i + j] == pblocks[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Block of size %" PetscInt_FMT " mismatches entry %" PetscInt_FMT, pblocks[i], pblocks[i + j]);
2843 }
2844 PetscCall(MatSetVariableBlockSizes(*J, nblocks, pblocks));
2845 }
2846 PetscCall(PetscFree(pblocks));
2847 }
2848 PetscCall(MatSetDM(*J, dm));
2849 PetscFunctionReturn(PETSC_SUCCESS);

--- 7695 unchanged lines hidden ---