plexfem.c (7307f5f486eb793bccb96258161f32863892bd17) plexfem.c (e91c04dfc8a52dee1965211bb1cc8e5bf775178f)
1#include "petscdm.h"
2#include "petscerror.h"
3#include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
4#include <petscsf.h>
5
6#include <petscblaslapack.h>
7#include <petsc/private/hashsetij.h>
8#include <petsc/private/petscfeimpl.h>

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

1621 const PetscInt debug = ((DM_Plex *)dm->data)->printL2;
1622 DM tdm;
1623 DMLabel depthLabel;
1624 PetscSection section;
1625 Vec localX, tv;
1626 PetscReal *localDiff;
1627 PetscInt dim, depth, dE, Nf, f, Nds, s;
1628 PetscBool transform;
1#include "petscdm.h"
2#include "petscerror.h"
3#include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
4#include <petscsf.h>
5
6#include <petscblaslapack.h>
7#include <petsc/private/hashsetij.h>
8#include <petsc/private/petscfeimpl.h>

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

1621 const PetscInt debug = ((DM_Plex *)dm->data)->printL2;
1622 DM tdm;
1623 DMLabel depthLabel;
1624 PetscSection section;
1625 Vec localX, tv;
1626 PetscReal *localDiff;
1627 PetscInt dim, depth, dE, Nf, f, Nds, s;
1628 PetscBool transform;
1629 PetscMPIInt Nfi;
1630
1631 PetscFunctionBegin;
1632 PetscCall(DMGetDimension(dm, &dim));
1633 PetscCall(DMGetCoordinateDim(dm, &dE));
1634 PetscCall(DMGetLocalSection(dm, &section));
1635 PetscCall(DMGetLocalVector(dm, &localX));
1636 PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
1637 PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));

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

1763 if (label) {
1764 PetscCall(ISRestoreIndices(pointIS, &points));
1765 PetscCall(ISDestroy(&pointIS));
1766 }
1767 PetscCall(ISRestoreIndices(fieldIS, &fields));
1768 PetscCall(PetscFree6(funcVal, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
1769 }
1770 PetscCall(DMRestoreLocalVector(dm, &localX));
1629
1630 PetscFunctionBegin;
1631 PetscCall(DMGetDimension(dm, &dim));
1632 PetscCall(DMGetCoordinateDim(dm, &dE));
1633 PetscCall(DMGetLocalSection(dm, &section));
1634 PetscCall(DMGetLocalVector(dm, &localX));
1635 PetscCall(DMGetBasisTransformDM_Internal(dm, &tdm));
1636 PetscCall(DMGetBasisTransformVec_Internal(dm, &tv));

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

1762 if (label) {
1763 PetscCall(ISRestoreIndices(pointIS, &points));
1764 PetscCall(ISDestroy(&pointIS));
1765 }
1766 PetscCall(ISRestoreIndices(fieldIS, &fields));
1767 PetscCall(PetscFree6(funcVal, interpolant, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
1768 }
1769 PetscCall(DMRestoreLocalVector(dm, &localX));
1771 PetscCall(PetscMPIIntCast(Nf, &Nfi));
1772 PetscCallMPI(MPIU_Allreduce(localDiff, diff, Nfi, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
1770 PetscCallMPI(MPIU_Allreduce(localDiff, diff, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
1773 PetscCall(PetscFree(localDiff));
1774 for (f = 0; f < Nf; ++f) diff[f] = PetscSqrtReal(diff[f]);
1775 PetscFunctionReturn(PETSC_SUCCESS);
1776}
1777
1778/*@C
1779 DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec.
1780

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

2522.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSNESComputeResidualFEM()`
2523@*/
2524PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
2525{
2526 PetscInt printFEM;
2527 PetscScalar *cintegral, *lintegral;
2528 PetscInt Nf, f, cellHeight, cStart, cEnd, cell;
2529 Vec locX;
1771 PetscCall(PetscFree(localDiff));
1772 for (f = 0; f < Nf; ++f) diff[f] = PetscSqrtReal(diff[f]);
1773 PetscFunctionReturn(PETSC_SUCCESS);
1774}
1775
1776/*@C
1777 DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec.
1778

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

2520.seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexSNESComputeResidualFEM()`
2521@*/
2522PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscScalar *integral, void *user)
2523{
2524 PetscInt printFEM;
2525 PetscScalar *cintegral, *lintegral;
2526 PetscInt Nf, f, cellHeight, cStart, cEnd, cell;
2527 Vec locX;
2530 PetscMPIInt Nfi;
2531
2532 PetscFunctionBegin;
2533 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2534 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
2535 PetscAssertPointer(integral, 3);
2536 PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
2537 PetscCall(DMPlexConvertPlex(dm, &dm, PETSC_TRUE));
2538 PetscCall(DMGetNumFields(dm, &Nf));

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

2550 printFEM = ((DM_Plex *)dm->data)->printFEM;
2551 /* Sum up values */
2552 for (cell = cStart; cell < cEnd; ++cell) {
2553 const PetscInt c = cell - cStart;
2554
2555 if (printFEM > 1) PetscCall(DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c * Nf]));
2556 for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c * Nf + f];
2557 }
2528
2529 PetscFunctionBegin;
2530 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2531 PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
2532 PetscAssertPointer(integral, 3);
2533 PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
2534 PetscCall(DMPlexConvertPlex(dm, &dm, PETSC_TRUE));
2535 PetscCall(DMGetNumFields(dm, &Nf));

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

2547 printFEM = ((DM_Plex *)dm->data)->printFEM;
2548 /* Sum up values */
2549 for (cell = cStart; cell < cEnd; ++cell) {
2550 const PetscInt c = cell - cStart;
2551
2552 if (printFEM > 1) PetscCall(DMPrintCellVector(cell, "Cell Integral", Nf, &cintegral[c * Nf]));
2553 for (f = 0; f < Nf; ++f) lintegral[f] += cintegral[c * Nf + f];
2554 }
2558 PetscCall(PetscMPIIntCast(Nf, &Nfi));
2559 PetscCallMPI(MPIU_Allreduce(lintegral, integral, Nfi, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
2555 PetscCallMPI(MPIU_Allreduce(lintegral, integral, Nf, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
2560 if (printFEM) {
2561 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Integral:"));
2562 for (f = 0; f < Nf; ++f) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " %g", (double)PetscRealPart(integral[f])));
2563 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "\n"));
2564 }
2565 PetscCall(PetscFree2(lintegral, cintegral));
2566 PetscCall(PetscLogEventEnd(DMPLEX_IntegralFEM, dm, 0, 0, 0));
2567 PetscCall(DMDestroy(&dm));

--- 4111 unchanged lines hidden --- 		2556 if (printFEM) {
2557 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Integral:"));
2558 for (f = 0; f < Nf; ++f) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " %g", (double)PetscRealPart(integral[f])));
2559 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "\n"));
2560 }
2561 PetscCall(PetscFree2(lintegral, cintegral));
2562 PetscCall(PetscLogEventEnd(DMPLEX_IntegralFEM, dm, 0, 0, 0));
2563 PetscCall(DMDestroy(&dm));

--- 4111 unchanged lines hidden ---