ohuih

miniapps/shifted/CMakeLists.txt

@@ -15,14 +15,24 @@ if (MFEM_USE_MPI)
         sbm_solver.cpp
         marking.cpp
         extrapolator.cpp
-	integ_algoim.cpp)
+        integ_algoim.cpp
+        shape_grad.cpp
+        mtop_solvers.cpp
+        mtop_filters.cpp
+	MMA.cpp)
+
+
   list(APPEND DIST_COMMON_HEADERS
         dist_solver.hpp
         sbm_solver.hpp
         sbm_aux.hpp
         marking.hpp
         extrapolator.hpp
-	integ_algoim.hpp)
+        integ_algoim.hpp
+        shape_grad.hpp
+        mtop_solvers.hpp
+        mtop_filters.hpp
+	MMA.hpp)
 
   convert_filenames_to_full_paths(DIST_COMMON_SOURCES)
   convert_filenames_to_full_paths(DIST_COMMON_HEADERS)
@@ -51,6 +61,40 @@ if (MFEM_USE_MPI)
     ${DIST_COMMON_FILES}
     LIBRARIES mfem mfem-common)
 
+  add_mfem_miniapp(shape_test
+    MAIN shape_grad_test.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(nodal_test
+      MAIN nodal_grad_test.cpp
+      ${DIST_COMMON_FILES}
+      LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(cfem_test
+      MAIN cfem_test.cpp
+      ${DIST_COMMON_FILES}
+      LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(volp_test
+      MAIN volp_test.cpp
+      ${DIST_COMMON_FILES}
+      LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(cobj_test
+      MAIN cut_obj_test.cpp
+      ${DIST_COMMON_FILES}
+      LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(exo2gmsh
+      MAIN exo2gmsh.cpp
+      ${DIST_COMMON_FILES}
+      LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(stress_bra
+      MAIN stress_bracket
+      LIBRARIES mfem mfem-common)
+
   if (MFEM_ENABLE_TESTING)
     add_test(NAME shifted_distance_np${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}

miniapps/shifted/MMA.cpp

@@ -0,0 +1,993 @@
+
+#include "MMA.hpp"
+#include <iostream>
+#include <math.h>
+
+/* -----------------------------------------------------------------------------
+Authors: Niels Aage
+ Copyright (C) 2013-2019,
+This MMA implementation is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+This Module is distributed in the hope that it will be useful,implementation
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+You should have received a copy of the GNU Lesser General Public
+License along with this Module; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+-------------------------------------------------------------------------- */
+
+MMA::MMA(PetscInt nn, PetscInt mm, PetscInt kk, Vec xo1t, Vec xo2t, Vec Ut, Vec Lt, PetscScalar* at, PetscScalar* ct,
+         PetscScalar* dt) {
+    n = nn;
+    m = mm;
+    k = kk;
+    if (k < 3) {
+        PetscPrintf(PETSC_COMM_WORLD, "NOT A LEGAL RESTART POINT (k<3): EXPECT BREAKDOWN\n");
+    }
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+    RobustAsymptotesType   = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    memcpy(a, at, mm * sizeof(PetscScalar));
+    memcpy(c, ct, mm * sizeof(PetscScalar));
+    memcpy(d, dt, mm * sizeof(PetscScalar));
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &L);
+    VecDuplicate(xo1t, &U);
+
+    VecDuplicate(xo1t, &alpha);
+    VecDuplicate(xo1t, &beta);
+
+    VecDuplicate(xo1t, &p0);
+    VecDuplicate(xo1t, &q0);
+    VecDuplicateVecs(xo1t, m, &pij);
+    VecDuplicateVecs(xo1t, m, &qij);
+
+    b = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &xo1);
+    VecDuplicate(xo1t, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+
+    // Now insert the values into xo1,xo2,U,L
+    PetscInt nloc;
+    VecGetLocalSize(xo1t, &nloc);
+
+    // input
+    PetscScalar *pxo1t, *pxo2t, *pUt, *pLt;
+    VecGetArray(xo1t, &pxo1t);
+    VecGetArray(xo2t, &pxo2t);
+    VecGetArray(Ut, &pUt);
+    VecGetArray(Lt, &pLt);
+
+    // internal
+    PetscScalar *pxo1, *pxo2, *pU, *pL;
+    VecGetArray(xo1, &pxo1);
+    VecGetArray(xo2, &pxo2);
+    VecGetArray(U, &pU);
+    VecGetArray(L, &pL);
+
+    // Copy data
+    memcpy(pxo1, pxo1t, nloc * sizeof(PetscScalar));
+    memcpy(pxo2, pxo2t, nloc * sizeof(PetscScalar));
+    memcpy(pU, pUt, nloc * sizeof(PetscScalar));
+    memcpy(pL, pLt, nloc * sizeof(PetscScalar));
+
+    // Restore arrays
+    VecRestoreArray(xo1t, &pxo1t);
+    VecRestoreArray(xo2t, &pxo2t);
+    VecRestoreArray(Ut, &pUt);
+    VecRestoreArray(Lt, &pLt);
+
+    VecRestoreArray(xo1, &pxo1);
+    VecRestoreArray(xo2, &pxo2);
+    VecRestoreArray(U, &pU);
+    VecRestoreArray(L, &pL);
+}
+
+MMA::MMA(PetscInt nn, PetscInt mm, PetscInt kk, Vec xo1t, Vec xo2t, Vec Ut, Vec Lt) {
+    n = nn;
+    m = mm;
+    k = kk;
+    if (k < 3) {
+        PetscPrintf(PETSC_COMM_WORLD, "NOT A LEGAL RESTART POINT (k<3): EXPECT BREAKDOWN\n");
+    }
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+    RobustAsymptotesType   = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    for (PetscInt i = 0; i < m; i++) {
+        a[i] = 0.0;
+        c[i] = 1000.0;
+        d[i] = 0.0;
+    }
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &L);
+    VecDuplicate(xo1t, &U);
+
+    VecDuplicate(xo1t, &alpha);
+    VecDuplicate(xo1t, &beta);
+
+    VecDuplicate(xo1t, &p0);
+    VecDuplicate(xo1t, &q0);
+    VecDuplicateVecs(xo1t, m, &pij);
+    VecDuplicateVecs(xo1t, m, &qij);
+
+    b = new PetscScalar[m];
+
+    VecDuplicate(xo1t, &xo1);
+    VecDuplicate(xo1t, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+
+    // Now insert the values into xo1,xo2,U,L
+    PetscInt nloc;
+    VecGetLocalSize(xo1t, &nloc);
+
+    // input
+    PetscScalar *pxo1t, *pxo2t, *pUt, *pLt;
+    VecGetArray(xo1t, &pxo1t);
+    VecGetArray(xo2t, &pxo2t);
+    VecGetArray(Ut, &pUt);
+    VecGetArray(Lt, &pLt);
+
+    // internal
+    PetscScalar *pxo1, *pxo2, *pU, *pL;
+    VecGetArray(xo1, &pxo1);
+    VecGetArray(xo2, &pxo2);
+    VecGetArray(U, &pU);
+    VecGetArray(L, &pL);
+
+    // Copy data
+    memcpy(pxo1, pxo1t, nloc * sizeof(PetscScalar));
+    memcpy(pxo2, pxo2t, nloc * sizeof(PetscScalar));
+    memcpy(pU, pUt, nloc * sizeof(PetscScalar));
+    memcpy(pL, pLt, nloc * sizeof(PetscScalar));
+
+    // Restore arrays
+    VecRestoreArray(xo1t, &pxo1t);
+    VecRestoreArray(xo2t, &pxo2t);
+    VecRestoreArray(Ut, &pUt);
+    VecRestoreArray(Lt, &pLt);
+
+    VecRestoreArray(xo1, &pxo1);
+    VecRestoreArray(xo2, &pxo2);
+    VecRestoreArray(U, &pU);
+    VecRestoreArray(L, &pL);
+}
+
+MMA::MMA(PetscInt nn, PetscInt mm, Vec x, PetscScalar* at, PetscScalar* ct, PetscScalar* dt) {
+    n = nn;
+    m = mm;
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+    RobustAsymptotesType   = 0;
+
+    k = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    memcpy(a, at, mm * sizeof(PetscScalar));
+    memcpy(c, ct, mm * sizeof(PetscScalar));
+    memcpy(d, dt, mm * sizeof(PetscScalar));
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(x, &L);
+    VecDuplicate(x, &U);
+
+    VecDuplicate(x, &alpha);
+    VecDuplicate(x, &beta);
+
+    VecDuplicate(x, &p0);
+    VecDuplicate(x, &q0);
+    VecDuplicateVecs(x, m, &pij);
+    VecDuplicateVecs(x, m, &qij);
+
+    b = new PetscScalar[m];
+
+    VecDuplicate(x, &xo1);
+    VecDuplicate(x, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+}
+
+MMA::MMA(PetscInt nn, PetscInt mm, Vec x) {
+    n = nn;
+    m = mm;
+
+    asyminit = 0.5;
+    asymdec  = 0.7;
+    asyminc  = 1.2;
+
+    NonLinConstraints      = PETSC_TRUE;
+    constraintModification = PETSC_FALSE;
+
+    RobustAsymptotesType = 0;
+
+    k = 0;
+
+    a = new PetscScalar[m];
+    c = new PetscScalar[m];
+    d = new PetscScalar[m];
+
+    for (PetscInt i = 0; i < m; i++) {
+        a[i] = 0.0;
+        c[i] = 1000.0;
+        d[i] = 0.0;
+    }
+
+    y   = new PetscScalar[m];
+    lam = new PetscScalar[m];
+
+    VecDuplicate(x, &L);
+    VecDuplicate(x, &U);
+
+    VecDuplicate(x, &alpha);
+    VecDuplicate(x, &beta);
+
+    VecDuplicate(x, &p0);
+    VecDuplicate(x, &q0);
+    VecDuplicateVecs(x, m, &pij);
+    VecDuplicateVecs(x, m, &qij);
+    b = new PetscScalar[m];
+
+    VecDuplicate(x, &xo1);
+    VecDuplicate(x, &xo2);
+
+    grad = new PetscScalar[m];
+    mu   = new PetscScalar[m];
+    s    = new PetscScalar[2 * m];
+
+    Hess = new PetscScalar[m * m];
+}
+
+MMA::~MMA() {
+
+    delete[] a;
+    delete[] b;
+    delete[] c;
+    delete[] d;
+    delete[] y;
+    delete[] lam;
+    VecDestroy(&L);
+    VecDestroy(&U);
+    VecDestroy(&alpha);
+    VecDestroy(&beta);
+    VecDestroy(&p0);
+    VecDestroy(&q0);
+    VecDestroyVecs(m, &pij);
+    VecDestroyVecs(m, &qij);
+    VecDestroy(&xo1);
+    VecDestroy(&xo2);
+    delete[] grad;
+    delete[] mu;
+    delete[] s;
+    delete[] Hess;
+}
+
+// restart method
+
+PetscErrorCode MMA::Restart(Vec xo1t, Vec xo2t, Vec Ut, Vec Lt) {
+
+    PetscErrorCode ierr = 0;
+
+    // Insert values into xo1t,xo2t,Ut,Lt
+    PetscInt nloc;
+    VecGetLocalSize(xo1t, &nloc);
+
+    // input
+    PetscScalar *pxo1t, *pxo2t, *pUt, *pLt;
+    VecGetArray(xo1t, &pxo1t);
+    VecGetArray(xo2t, &pxo2t);
+    VecGetArray(Ut, &pUt);
+    VecGetArray(Lt, &pLt);
+
+    // internal
+    PetscScalar *pxo1, *pxo2, *pU, *pL;
+    VecGetArray(xo1, &pxo1);
+    VecGetArray(xo2, &pxo2);
+    VecGetArray(U, &pU);
+    VecGetArray(L, &pL);
+
+    // Copy data
+    memcpy(pxo1t, pxo1, nloc * sizeof(PetscScalar));
+    memcpy(pxo2t, pxo2, nloc * sizeof(PetscScalar));
+    memcpy(pUt, pU, nloc * sizeof(PetscScalar));
+    memcpy(pLt, pL, nloc * sizeof(PetscScalar));
+
+    // Restore arrays
+    VecRestoreArray(xo1t, &pxo1t);
+    VecRestoreArray(xo2t, &pxo2t);
+    VecRestoreArray(Ut, &pUt);
+    VecRestoreArray(Lt, &pLt);
+
+    VecRestoreArray(xo1, &pxo1);
+    VecRestoreArray(xo2, &pxo2);
+    VecRestoreArray(U, &pU);
+    VecRestoreArray(L, &pL);
+
+    return (ierr);
+}
+
+// Set the aggresivity of the moving asymptotes
+PetscErrorCode MMA::SetAsymptotes(PetscScalar init, PetscScalar decrease, PetscScalar increase) {
+    PetscErrorCode ierr = 0;
+
+    // asymptotes initialization and increase/decrease
+    asyminit = init;
+    asymdec  = decrease;
+    asyminc  = increase;
+    return ierr;
+}
+
+PetscErrorCode MMA::SetRobustAsymptotesType(PetscInt val) {
+
+    PetscErrorCode ierr = 0;
+
+    RobustAsymptotesType = val;
+
+    if (RobustAsymptotesType == 0 || RobustAsymptotesType == 1) {
+    } else {
+        RobustAsymptotesType = 0;
+        PetscPrintf(PETSC_COMM_WORLD, "ERROR in MMA.cc/h: RobustAsymptotesType cannot be set to: %d \n", val);
+    }
+    return ierr;
+}
+
+PetscErrorCode MMA::SetOuterMovelimit(PetscScalar Xmin, PetscScalar Xmax, PetscScalar movlim, Vec x, Vec xmin,
+                                      Vec xmax) {
+
+    PetscErrorCode ierr = 0;
+
+    PetscScalar *xv, *xmiv, *xmav;
+    PetscInt     nloc;
+    VecGetLocalSize(x, &nloc);
+    VecGetArray(x, &xv);
+    VecGetArray(xmin, &xmiv);
+    VecGetArray(xmax, &xmav);
+    for (PetscInt i = 0; i < nloc; i++) {
+        xmav[i] = Min(Xmax, xv[i] + movlim);
+        xmiv[i] = Max(Xmin, xv[i] - movlim);
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArray(xmin, &xmiv);
+    VecRestoreArray(xmax, &xmav);
+    return ierr;
+}
+
+PetscScalar MMA::DesignChange(Vec x, Vec xold) {
+
+    PetscScalar *xv, *xo;
+    PetscInt     nloc;
+    VecGetLocalSize(x, &nloc);
+    VecGetArray(x, &xv);
+    VecGetArray(xold, &xo);
+    PetscScalar ch = 0.0;
+    for (PetscInt i = 0; i < nloc; i++) {
+        ch    = PetscMax(ch, PetscAbsReal(xv[i] - xo[i]));
+        xo[i] = xv[i];
+    }
+    PetscScalar tmp;
+    MPI_Allreduce(&ch, &tmp, 1, MPIU_SCALAR, MPI_MAX, PETSC_COMM_WORLD);
+    ch = tmp;
+    VecRestoreArray(x, &xv);
+    VecRestoreArray(xold, &xo);
+
+    return (ch);
+}
+
+PetscErrorCode MMA::KKTresidual(Vec x, Vec dfdx, PetscScalar* fx, Vec* dgdx, Vec xmin, Vec xmax, PetscScalar* norm2,
+                                PetscScalar* normInf) {
+    PetscErrorCode ierr = 0;
+
+    if (!NonLinConstraints) {
+        PetscErrorPrintf("MMA->KKTresidual called WITH constraints but object was "
+                         "allocated WITHOUT !\n");
+        return -1;
+    }
+
+    PetscScalar *xp, *xminp, *xmaxp, *df0dxp, **dfdxp;
+    PetscInt     locsiz;
+    VecGetLocalSize(x, &locsiz);
+
+    VecGetArray(x, &xp);
+    VecGetArray(xmin, &xminp);
+    VecGetArray(xmax, &xmaxp);
+    VecGetArray(dfdx, &df0dxp);
+    VecGetArrays(dgdx, m, &dfdxp);
+
+    PetscScalar resi, ri, mu_min, mu_max;
+
+    norm2[0]   = 0;
+    normInf[0] = 0;
+    for (PetscInt i = 0; i < locsiz; i++) {
+        ri = df0dxp[i];
+        for (PetscInt j = 0; j < m; j++) {
+            ri += lam[j] * dfdxp[j][i];
+        }
+        mu_min = 0.0;
+        if (xp[i] < xminp[i] + 1.0e-5 && ri > 0.0) {
+            mu_min = ri;
+        }
+        mu_max = 0.0;
+        if (xp[i] > xmaxp[i] - 1.0e-5 && ri < 0.0) {
+            mu_max = -ri;
+        }
+        ri += -mu_min + mu_max;
+        norm2[0] += pow(ri, 2.0);
+        normInf[0] = Max(Abs(ri), normInf[0]);
+        resi       = mu_min * (xp[i] - xminp[i]);
+        norm2[0] += pow(resi, 2.0);
+        normInf[0] = Max(Abs(resi), normInf[0]);
+        resi       = mu_max * (xmaxp[i] - xp[i]);
+        norm2[0] += pow(resi, 2.0);
+        normInf[0] = Max(Abs(resi), normInf[0]);
+    }
+
+    VecRestoreArray(x, &xp);
+    VecRestoreArray(xmin, &xminp);
+    VecRestoreArray(xmax, &xmaxp);
+    VecRestoreArray(dfdx, &df0dxp);
+    VecRestoreArrays(dgdx, m, &dfdxp);
+    PetscScalar n2tmp = norm2[0];
+    PetscScalar nItmp = normInf[0];
+    norm2[0]          = 0.0;
+    normInf[0]        = 0.0;
+    MPI_Allreduce(&n2tmp, norm2, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+    MPI_Allreduce(&nItmp, normInf, 1, MPIU_SCALAR, MPI_MAX, PETSC_COMM_WORLD);
+    ri = 0.0;
+    for (PetscInt j = 0; j < m; j++) {
+        ri += lam[j] * (a[j] * z + y[j] - fx[j]);
+    }
+    norm2[0] += pow(ri, 2.0);
+    normInf[0] = Max(Abs(ri), normInf[0]);
+    norm2[0]   = sqrt(norm2[0]);
+
+    return ierr;
+}
+
+// Set and solve a subproblem: return new xval
+PetscErrorCode MMA::Update(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax) {
+    PetscErrorCode ierr = 0;
+
+    if (!NonLinConstraints) {
+        PetscErrorPrintf("MMA->Update called WITH constraints but object was "
+                         "allocated WITHOUT !\n");
+        return -1;
+    }
+
+    // Generate the subproblem
+    GenSub(xval, dfdx, gx, dgdx, xmin, xmax);
+
+    // Update xolds
+    VecCopy(xo1, xo2);
+    VecCopy(xval, xo1);
+
+    // Solve the dual with an interior point method
+    SolveDIP(xval);
+    return ierr;
+}
+
+// PRIVATE METHODS
+
+PetscErrorCode MMA::GenSub(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax) {
+    PetscErrorCode ierr = 0;
+
+    PetscScalar gamma, helpvar;
+
+    k++;
+    PetscInt nloc;
+    VecGetLocalSize(xval, &nloc);
+    PetscScalar *xv, *Lv, *Uv, *x1v, *x2v, *xminv, *xmaxv;
+    PetscScalar *alf, *bet, *dfdxv, *p0v, *q0v, **dgdxv, **pijv, **qijv;
+    if (k < 3) {
+        VecAXPBYPCZ(L, (PetscScalar)1.0, -asyminit, (PetscScalar)0.0, xval, xmax);
+        VecAXPY(L, asyminit, xmin);
+        VecAXPBYPCZ(U, (PetscScalar)1.0, +asyminit, (PetscScalar)0.0, xval, xmax);
+        VecAXPY(U, -asyminit, xmin);
+    }
+    VecGetArray(xval, &xv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    VecGetArray(xo1, &x1v);
+    VecGetArray(xo2, &x2v);
+    VecGetArray(xmin, &xminv);
+    VecGetArray(xmax, &xmaxv);
+
+    VecGetArray(alpha, &alf);
+    VecGetArray(beta, &bet);
+    VecGetArray(dfdx, &dfdxv);
+    VecGetArray(p0, &p0v);
+    VecGetArray(q0, &q0v);
+
+    VecGetArrays(dgdx, m, &dgdxv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    if (k > 2) {
+        for (PetscInt i = 0; i < nloc; i++) {
+            helpvar = (xv[i] - x1v[i]) * (x1v[i] - x2v[i]);
+            if (helpvar < 0.0) {
+                gamma = asymdec;
+            } else if (helpvar > 0.0) {
+                gamma = asyminc;
+            } else {
+                gamma = 1.0;
+            }
+            Lv[i] = xv[i] - gamma * (x1v[i] - Lv[i]);
+            Uv[i] = xv[i] + gamma * (Uv[i] - x1v[i]);
+            PetscScalar xmi, xma;
+            xmi = Max(1.0e-5, xmaxv[i] - xminv[i]);
+            if (RobustAsymptotesType == 0) {
+                Lv[i] = Max(Lv[i], xv[i] - 10.0 * xmi);
+                Lv[i] = Min(Lv[i], xv[i] - 0.01 * xmi);
+                Uv[i] = Max(Uv[i], xv[i] + 0.01 * xmi);
+                Uv[i] = Min(Uv[i], xv[i] + 10.0 * xmi);
+            } else if (RobustAsymptotesType == 1) {
+                Lv[i] = Max(Lv[i], xv[i] - 100.0 * xmi);
+                Lv[i] = Min(Lv[i], xv[i] - 1.0e-4 * xmi);
+                Uv[i] = Max(Uv[i], xv[i] + 1.0e-4 * xmi);
+                Uv[i] = Min(Uv[i], xv[i] + 100.0 * xmi);
+                xmi   = xminv[i] - 1.0e-5;
+                xma   = xmaxv[i] + 1.0e-5;
+                if (xv[i] < xmi) {
+                    Lv[i] = xv[i] - (xma - xv[i]) / 0.9;
+                    Uv[i] = xv[i] + (xma - xv[i]) / 0.9;
+                }
+                if (xv[i] > xma) {
+                    Lv[i] = xv[i] - (xv[i] - xmi) / 0.9;
+                    Uv[i] = xv[i] + (xv[i] - xmi) / 0.9;
+                }
+            }
+        }
+    }
+    PetscScalar dfdxp, dfdxm;
+    PetscScalar feps = 1.0e-6;
+    for (PetscInt i = 0; i < nloc; i++) {
+        alf[i] = Max(xminv[i], 0.9 * Lv[i] + 0.1 * xv[i]);
+        bet[i] = Min(xmaxv[i], 0.9 * Uv[i] + 0.1 * xv[i]);
+        dfdxp  = Max(0.0, dfdxv[i]);
+        dfdxm  = Max(0.0, -1.0 * dfdxv[i]);
+        p0v[i] = pow(Uv[i] - xv[i], 2.0) * (dfdxp + 0.001 * Abs(dfdxv[i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+        q0v[i] = pow(xv[i] - Lv[i], 2.0) * (dfdxm + 0.001 * Abs(dfdxv[i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+        for (PetscInt j = 0; j < m; j++) {
+            dfdxp = Max(0.0, dgdxv[j][i]);
+            dfdxm = Max(0.0, -1.0 * dgdxv[j][i]);
+            if (constraintModification) {
+                pijv[j][i] =
+                    pow(Uv[i] - xv[i], 2.0) * (dfdxp + 0.001 * Abs(dgdxv[j][i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+                qijv[j][i] =
+                    pow(xv[i] - Lv[i], 2.0) * (dfdxm + 0.001 * Abs(dgdxv[j][i]) + 0.5 * feps / (Uv[i] - Lv[i]));
+            } else {
+                pijv[j][i] = pow(Uv[i] - xv[i], 2.0) * (dfdxp);
+                qijv[j][i] = pow(xv[i] - Lv[i], 2.0) * (dfdxm);
+            }
+        }
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        b[j] = 0.0;
+        for (PetscInt i = 0; i < nloc; i++) {
+            b[j] += pijv[j][i] / (Uv[i] - xv[i]) + qijv[j][i] / (xv[i] - Lv[i]);
+        }
+    }
+    {
+        PetscScalar* tmp = new PetscScalar[m];
+        for (PetscInt i = 0; i < m; i++) {
+            tmp[i] = 0.0;
+        }
+        MPI_Allreduce(b, tmp, m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(b, tmp, sizeof(PetscScalar) * m);
+        delete[] tmp;
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        b[j] += -gx[j];
+    }
+    VecRestoreArray(xval, &xv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    VecRestoreArray(xo1, &x1v);
+    VecRestoreArray(xo2, &x2v);
+    VecRestoreArray(xmin, &xminv);
+    VecRestoreArray(xmax, &xmaxv);
+
+    VecRestoreArray(alpha, &alf);
+    VecRestoreArray(beta, &bet);
+    VecRestoreArray(dfdx, &dfdxv);
+
+    VecRestoreArrays(dgdx, m, &dgdxv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    return ierr;
+}
+
+PetscErrorCode MMA::SolveDIP(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    for (PetscInt j = 0; j < m; j++) {
+        lam[j] = c[j] / 2.0;
+        mu[j]  = 1.0;
+    }
+    PetscScalar tol  = 1.0e-9 * sqrt(m + n);
+    PetscScalar epsi = 1.0;
+    PetscScalar err  = 1.0;
+    PetscInt    loop;
+    while (epsi > tol) {
+
+        loop = 0;
+        while (err > 0.9 * epsi && loop < 100) {
+            loop++;
+            XYZofLAMBDA(x);
+            DualGrad(x);
+            for (PetscInt j = 0; j < m; j++) {
+                grad[j] = -1.0 * grad[j] - epsi / lam[j];
+            }
+            DualHess(x);
+            Factorize(Hess, m);
+            Solve(Hess, grad, m);
+            for (PetscInt j = 0; j < m; j++) {
+                s[j] = grad[j];
+            }
+            for (PetscInt i = 0; i < m; i++) {
+                s[m + i] = -mu[i] + epsi / lam[i] - s[i] * mu[i] / lam[i];
+            }
+            DualLineSearch();
+            XYZofLAMBDA(x);
+            err = DualResidual(x, epsi);
+        }
+        epsi = epsi * 0.1;
+    }
+    return ierr;
+}
+
+PetscErrorCode MMA::XYZofLAMBDA(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar *xv, **pijv, **qijv, *p0v, *q0v, *alf, *bet, *Lv, *Uv;
+    VecGetArray(x, &xv);
+    VecGetArray(p0, &p0v);
+    VecGetArray(q0, &q0v);
+    VecGetArray(alpha, &alf);
+    VecGetArray(beta, &bet);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    PetscScalar lamai = 0.0;
+    for (PetscInt i = 0; i < m; i++) {
+        if (lam[i] < 0.0) {
+            lam[i] = 0;
+        }
+        y[i] = Max(0.0, lam[i] - c[i]);
+        lamai += lam[i] * a[i];
+    }
+    z = Max(0.0, 10.0 * (lamai - 1.0)); // SINCE a0 = 1.0
+    PetscScalar pjlam, qjlam;
+    for (PetscInt i = 0; i < nloc; i++) {
+        pjlam = p0v[i];
+        qjlam = q0v[i];
+        for (PetscInt j = 0; j < m; j++) {
+            pjlam += pijv[j][i] * lam[j];
+            qjlam += qijv[j][i] * lam[j];
+        }
+        xv[i] = (sqrt(pjlam) * Lv[i] + sqrt(qjlam) * Uv[i]) / (sqrt(pjlam) + sqrt(qjlam));
+        if (xv[i] < alf[i]) {
+            xv[i] = alf[i];
+        }
+        if (xv[i] > bet[i]) {
+            xv[i] = bet[i];
+        }
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(p0, &p0v);
+    VecRestoreArray(q0, &q0v);
+    VecRestoreArray(alpha, &alf);
+    VecRestoreArray(beta, &bet);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    return ierr;
+}
+
+PetscErrorCode MMA::DualGrad(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar *xv, *Lv, *Uv, **pijv, **qijv;
+    VecGetArray(x, &xv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    for (PetscInt j = 0; j < m; j++) {
+        grad[j] = 0.0;
+        for (PetscInt i = 0; i < nloc; i++) {
+            grad[j] += pijv[j][i] / (Uv[i] - xv[i]) + qijv[j][i] / (xv[i] - Lv[i]);
+        }
+    }
+    {
+        PetscScalar* tmp = new PetscScalar[m];
+        for (PetscInt i = 0; i < m; i++) {
+            tmp[i] = 0.0;
+        }
+        MPI_Allreduce(grad, tmp, m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(grad, tmp, sizeof(PetscScalar) * m);
+        delete[] tmp;
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        grad[j] += -b[j] - a[j] * z - y[j];
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    return ierr;
+}
+
+PetscErrorCode MMA::DualHess(Vec x) {
+    PetscErrorCode ierr = 0;
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar *xv, *Lv, *Uv, **pijv, **qijv, *alf, *bet, *p0v, *q0v;
+    VecGetArray(x, &xv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    VecGetArray(alpha, &alf);
+    VecGetArray(beta, &bet);
+    VecGetArray(p0, &p0v);
+    VecGetArray(q0, &q0v);
+    PetscScalar* df2 = new PetscScalar[nloc];
+    PetscScalar* PQ  = new PetscScalar[nloc * m];
+    PetscScalar  pjlam, qjlam;
+    for (PetscInt i = 0; i < nloc; i++) {
+        pjlam = p0v[i];
+        qjlam = q0v[i];
+        for (PetscInt j = 0; j < m; j++) {
+            pjlam += pijv[j][i] * lam[j];
+            qjlam += qijv[j][i] * lam[j];
+            PQ[i * m + j] = pijv[j][i] / pow(Uv[i] - xv[i], 2.0) - qijv[j][i] / pow(xv[i] - Lv[i], 2.0);
+        }
+        df2[i]         = -1.0 / (2.0 * pjlam / pow(Uv[i] - xv[i], 3.0) + 2.0 * qjlam / pow(xv[i] - Lv[i], 3.0));
+        PetscScalar xp = (sqrt(pjlam) * Lv[i] + sqrt(qjlam) * Uv[i]) / (sqrt(pjlam) + sqrt(qjlam));
+        if (xp < alf[i]) {
+            df2[i] = 0.0;
+        }
+        if (xp > bet[i]) {
+            df2[i] = 0.0;
+        }
+    }
+    PetscScalar* tmp = new PetscScalar[n * m];
+    for (PetscInt j = 0; j < m; j++) {
+        for (PetscInt i = 0; i < nloc; i++) {
+            tmp[j * nloc + i] = 0.0;
+            tmp[j * nloc + i] += PQ[i * m + j] * df2[i];
+        }
+    }
+    for (PetscInt i = 0; i < m; i++) {
+        for (PetscInt j = 0; j < m; j++) {
+            Hess[i * m + j] = 0.0;
+            for (PetscInt k = 0; k < nloc; k++) {
+                Hess[i * m + j] += tmp[i * nloc + k] * PQ[k * m + j];
+            }
+        }
+    }
+    {
+        PetscScalar* tmpp = new PetscScalar[m * m];
+        for (PetscInt i = 0; i < m * m; i++) {
+            tmpp[i] = Hess[i];
+        }
+        MPI_Allreduce(Hess, tmpp, m * m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(Hess, tmpp, sizeof(PetscScalar) * m * m);
+        delete[] tmpp;
+    }
+    PetscScalar lamai = 0.0;
+    for (PetscInt j = 0; j < m; j++) {
+        if (lam[j] < 0.0) {
+            lam[j] = 0.0;
+        }
+        lamai += lam[j] * a[j];
+        if (lam[j] > c[j]) {
+            Hess[j * m + j] += -1.0;
+        }
+        Hess[j * m + j] += -mu[j] / lam[j];
+    }
+    if (lamai > 0.0) {
+        for (PetscInt j = 0; j < m; j++) {
+            for (PetscInt k = 0; k < m; k++) {
+                Hess[j * m + k] += -10.0 * a[j] * a[k];
+            }
+        }
+    }
+    PetscScalar HessTrace = 0.0;
+    for (PetscInt i = 0; i < m; i++) {
+        HessTrace += Hess[i * m + i];
+    }
+    PetscScalar HessCorr = 1e-4 * HessTrace / m;
+    if (-1.0 * HessCorr < 1.0e-7) {
+        HessCorr = -1.0e-7;
+    }
+    for (PetscInt i = 0; i < m; i++) {
+        Hess[i * m + i] += HessCorr;
+    }
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    VecRestoreArray(q0, &q0v);
+    VecRestoreArray(p0, &p0v);
+    VecRestoreArray(alpha, &alf);
+    VecRestoreArray(beta, &bet);
+    delete[] df2;
+    delete[] PQ;
+    delete[] tmp;
+    return ierr;
+}
+
+PetscErrorCode MMA::DualLineSearch() {
+    PetscErrorCode ierr = 0;
+
+    PetscScalar theta = 1.005;
+    for (PetscInt i = 0; i < m; i++) {
+        if (theta < -1.01 * s[i] / lam[i]) {
+            theta = -1.01 * s[i] / lam[i];
+        }
+        if (theta < -1.01 * s[i + m] / mu[i]) {
+            theta = -1.01 * s[i + m] / mu[i];
+        }
+    }
+    theta = 1.0 / theta;
+    for (PetscInt i = 0; i < m; i++) {
+        lam[i] = lam[i] + theta * s[i];
+        mu[i]  = mu[i] + theta * s[i + m];
+    }
+    return ierr;
+}
+
+PetscScalar MMA::DualResidual(Vec x, PetscScalar epsi) {
+
+    PetscInt nloc;
+    VecGetLocalSize(x, &nloc);
+    PetscScalar* res = new PetscScalar[2 * m];
+    PetscScalar *xv, *Lv, *Uv, **pijv, **qijv;
+    VecGetArray(x, &xv);
+    VecGetArrays(pij, m, &pijv);
+    VecGetArrays(qij, m, &qijv);
+    VecGetArray(L, &Lv);
+    VecGetArray(U, &Uv);
+    for (PetscInt j = 0; j < m; j++) {
+        res[j]     = 0.0;
+        res[j + m] = 0.0;
+        for (PetscInt i = 0; i < nloc; i++) {
+            res[j] += pijv[j][i] / (Uv[i] - xv[i]) + qijv[j][i] / (xv[i] - Lv[i]);
+        }
+    }
+    {
+        PetscScalar* tmp = new PetscScalar[2 * m];
+        for (PetscInt i = 0; i < 2 * m; i++) {
+            tmp[i] = 0.0;
+        }
+        MPI_Allreduce(res, tmp, 2 * m, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
+        memcpy(res, tmp, sizeof(PetscScalar) * 2 * m);
+        delete[] tmp;
+    }
+    for (PetscInt j = 0; j < m; j++) {
+        res[j] += -b[j] - a[j] * z - y[j] + mu[j];
+        res[j + m] += mu[j] * lam[j] - epsi;
+    }
+    PetscScalar nrI = 0.0;
+    for (PetscInt i = 0; i < 2 * m; i++) {
+        if (nrI < Abs(res[i])) {
+            nrI = Abs(res[i]);
+        }
+    }
+    delete[] res;
+    VecRestoreArray(x, &xv);
+    VecRestoreArrays(pij, m, &pijv);
+    VecRestoreArrays(qij, m, &qijv);
+    VecRestoreArray(L, &Lv);
+    VecRestoreArray(U, &Uv);
+    return nrI;
+}
+
+PetscErrorCode MMA::Factorize(PetscScalar* K, PetscInt nn) {
+    PetscErrorCode ierr = 0;
+
+    for (PetscInt ss = 0; ss < nn - 1; ss++) {
+        for (PetscInt i = ss + 1; i < nn; i++) {
+            K[i * nn + ss] = K[i * nn + ss] / K[ss * nn + ss];
+            for (PetscInt j = ss + 1; j < nn; j++) {
+                K[i * nn + j] = K[i * nn + j] - K[i * nn + ss] * K[ss * nn + j];
+            }
+        }
+    }
+    return ierr;
+}
+
+PetscErrorCode MMA::Solve(PetscScalar* K, PetscScalar* x, PetscInt nn) {
+    PetscErrorCode ierr = 0;
+
+    for (PetscInt i = 1; i < nn; i++) {
+        PetscScalar a = 0.0;
+        for (PetscInt j = 0; j < i; j++) {
+            a = a - K[i * nn + j] * x[j];
+        }
+        x[i] = x[i] + a;
+    }
+    x[nn - 1] = x[nn - 1] / K[(nn - 1) * nn + (nn - 1)];
+    for (PetscInt i = nn - 2; i >= 0; i--) {
+        PetscScalar a = x[i];
+        for (PetscInt j = i + 1; j < nn; j++) {
+            a = a - K[i * nn + j] * x[j];
+        }
+        x[i] = a / K[i * nn + i];
+    }
+    return ierr;
+}
+
+PetscScalar MMA::Min(PetscScalar d1, PetscScalar d2) { return d1 < d2 ? d1 : d2; }
+
+PetscScalar MMA::Max(PetscScalar d1, PetscScalar d2) { return d1 > d2 ? d1 : d2; }
+
+PetscInt MMA::Min(PetscInt d1, PetscInt d2) { return d1 < d2 ? d1 : d2; }
+
+PetscInt MMA::Max(PetscInt d1, PetscInt d2) { return d1 > d2 ? d1 : d2; }
+
+PetscScalar MMA::Abs(PetscScalar d1) { return d1 > 0 ? d1 : -1.0 * d1; }
+
+

miniapps/shifted/MMA.hpp

@@ -0,0 +1,240 @@
+#ifndef __MMA__HPP
+#define __MMA__HPP
+
+#include "mfem.hpp"
+#include <petsc.h>
+
+/* -----------------------------------------------------------------------------
+Authors: Niels Aage
+ Copyright (C) 2013-2020,
+This MMA implementation is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+This Module is distributed in the hope that it will be useful,implementation
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+You should have received a copy of the GNU Lesser General Public
+License along with this Module; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+-------------------------------------------------------------------------- */
+class MMA {
+  public:
+    // Construct using defaults subproblem penalization
+    MMA(PetscInt n, PetscInt m, Vec x);
+    // User defined subproblem penalization
+    MMA(PetscInt n, PetscInt m, Vec x, PetscScalar* a, PetscScalar* c, PetscScalar* d);
+    // Initialize with restart from itr
+    MMA(PetscInt n, PetscInt m, PetscInt itr, Vec xo1, Vec xo2, Vec U, Vec L);
+    // Initialize with restart and specify subproblem parameters
+    MMA(PetscInt n, PetscInt m, PetscInt itr, Vec xo1, Vec xo2, Vec U, Vec L, PetscScalar* a, PetscScalar* c,
+        PetscScalar* d);
+    // Destructor
+    ~MMA();
+
+    // Set and solve a subproblem: return new xval
+    PetscErrorCode Update(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax);
+
+    // Return necessary data for possible restart
+    PetscErrorCode Restart(Vec xo1, Vec xo2, Vec U, Vec L);
+
+    // Set the aggresivity of the moving asymptotes
+    PetscErrorCode SetAsymptotes(PetscScalar init, PetscScalar decrease, PetscScalar increase);
+
+    // do/don't add convexity approx to constraints: default=false
+    PetscErrorCode ConstraintModification(PetscBool conMod) {
+        constraintModification = conMod;
+        return 0;
+    };
+
+    // val=0: default, val=1: increase robustness, i.e
+    // control the spacing between L < alp < x < beta < U,
+    PetscErrorCode SetRobustAsymptotesType(PetscInt val);
+
+    // Sets outer movelimits on all primal design variables
+    // This is often requires to prevent the solver from oscilating
+    PetscErrorCode SetOuterMovelimit(PetscScalar Xmin, PetscScalar Xmax, PetscScalar movelim, Vec x, Vec xmin,
+                                     Vec xmax);
+
+    // Return KKT residual norms (norm2 and normInf)
+    PetscErrorCode KKTresidual(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax, PetscScalar* norm2,
+                               PetscScalar* normInf);
+
+    // Inf norm on diff between two vectors: SHOULD NOT BE HERE - USE BASIC
+    // PETSc!!!!!
+    PetscScalar DesignChange(Vec x, Vec xold);
+
+  private:
+    // Set up the MMA subproblem based on old x's and xval
+    PetscErrorCode GenSub(Vec xval, Vec dfdx, PetscScalar* gx, Vec* dgdx, Vec xmin, Vec xmax);
+
+    // Interior point solver for the subproblem
+    PetscErrorCode SolveDIP(Vec xval);
+
+    // Compute primal vars based on dual solution
+    PetscErrorCode XYZofLAMBDA(Vec x);
+
+    // Dual gradient
+    PetscErrorCode DualGrad(Vec x);
+
+    // Dual Hessian
+    PetscErrorCode DualHess(Vec x);
+
+    // Dual line search
+    PetscErrorCode DualLineSearch();
+
+    // Dual residual
+    PetscScalar DualResidual(Vec x, PetscScalar epsi);
+
+    // Problem size and iteration counter
+    PetscInt n, m, k;
+
+    // "speed-control" for the asymptotes
+    PetscScalar asyminit, asymdec, asyminc;
+
+    // do/don't add convexity constraint approximation in subproblem
+    PetscBool constraintModification; // default = FALSE
+
+    // Bool specifying if non lin constraints are included or not
+    PetscBool NonLinConstraints;
+
+    // 0: (default) span between alp L x U beta,
+    // 1: increase the span for further robustness
+    PetscInt RobustAsymptotesType;
+
+    // Local vectors: penalty numbers for subproblem
+    PetscScalar *a, *c, *d;
+
+    // Local vectors: elastic variables
+    PetscScalar* y;
+    PetscScalar  z;
+
+    // Local vectors: Lagrange multipliers:
+    PetscScalar *lam, *mu, *s;
+
+    // Global: Asymptotes, bounds, objective approx., constraint approx.
+    Vec L, U, alpha, beta, p0, q0, *pij, *qij;
+
+    // Local: subproblem constant terms, dual gradient, dual hessian
+    PetscScalar *b, *grad, *Hess;
+
+    // Global: Old design variables
+    Vec xo1, xo2;
+
+    // Math helpers
+    PetscErrorCode Factorize(PetscScalar* K, PetscInt nn);
+    PetscErrorCode Solve(PetscScalar* K, PetscScalar* x, PetscInt nn);
+    PetscScalar    Min(PetscScalar d1, PetscScalar d2);
+    PetscScalar    Max(PetscScalar d1, PetscScalar d2);
+    PetscInt       Min(PetscInt d1, PetscInt d2);
+    PetscInt       Max(PetscInt d1, PetscInt d2);
+    PetscScalar    Abs(PetscScalar d1);
+};
+
+
+namespace mfem {
+
+class NativeMMA
+{
+public:
+    // User defined subproblem penalization
+    NativeMMA(MPI_Comm comm_, int m, mfem::Vector& x, double* a, double* c, double* d){
+
+        comm=comm_;
+        num_con=m;
+        Vec pv;
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&pv);
+
+        PetscScalar* ap=new PetscScalar[m];
+        PetscScalar* cp=new PetscScalar[m];
+        PetscScalar* dp=new PetscScalar[m];
+        for(int i=0;i<m;i++)
+        {
+            ap[i]=a[i];
+            cp[i]=c[i];
+            dp[i]=c[i];
+        }
+
+        PetscInt nn;
+        VecGetSize(pv,&nn);
+
+        mma=new MMA(nn, m, pv, ap, cp, dp);
+
+        delete [] ap;
+        delete [] cp;
+        delete [] dp;
+
+        VecDestroy(&pv);
+
+        //allocate the native PETSc objects necessary for the subproblems
+
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&xval);
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&dfdx);
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&xmin);
+        VecCreateMPI(comm,x.Size(),PETSC_DETERMINE,&xmax);
+
+        VecDuplicateVecs(xval,num_con, &dgdx);
+    }
+
+    ~NativeMMA()
+    {
+        delete mma;
+
+        VecDestroy(&xval);
+        VecDestroy(&dfdx);
+        VecDestroy(&xmin);
+        VecDestroy(&xmax);
+        VecDestroyVecs(num_con,&dgdx);
+    }
+
+    // Set and solve a subproblem: return new xval
+    void Update(Vector& xval_, Vector& dfdx_, double* gx_, Vector* dgdx_, Vector& xmin_, Vector& xmax_){
+        //copy data
+        double* data;
+        VecGetArray(xval,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=xval_[i];}
+        VecRestoreArray(xval,&data);
+        //dfdx
+        VecGetArray(dfdx,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=dfdx_[i];}
+        VecRestoreArray(dfdx,&data);
+        //dgdx
+        for(int j=0;j<num_con;j++){
+            VecGetArray(dgdx[j],&data);
+            for(int i=0;i<xval_.Size();i++){data[i]=(dgdx_[j])[i];}
+            VecRestoreArray(dgdx[j],&data);
+        }
+        //xmin
+        VecGetArray(xmin,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=xmin_[i];}
+        VecRestoreArray(xmin,&data);
+        //xmax
+        VecGetArray(xmax,&data);
+        for(int i=0;i<xval_.Size();i++){data[i]=xmax_[i];}
+        VecRestoreArray(xmax,&data);
+
+        mma->Update(xval,dfdx,gx_,dgdx,xmin,xmax);
+
+        VecGetArray(xval,&data);
+        for(int i=0;i<xval_.Size();i++){xval_[i]=data[i];}
+        VecRestoreArray(xval,&data);
+    }
+private:
+    MMA* mma;
+    int num_con;
+    MPI_Comm comm;
+
+    Vec xval;
+    Vec dfdx;
+    Vec* dgdx;
+    Vec xmin;
+    Vec xmax;
+
+
+};
+
+}
+
+
+#endif

miniapps/shifted/cfem_test.cpp

@@ -0,0 +1,676 @@
+#include <fstream>
+#include <iostream>
+
+#include "mfem.hpp"
+#include "marking.hpp"
+#include "mtop_solvers.hpp"
+#include "mtop_filters.hpp"
+
+
+using namespace mfem;
+using namespace std;
+
+class GyroidCoeff:public Coefficient
+{
+public:
+    GyroidCoeff(double cell_size=1.0){
+        ll=cell_size;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+        double x = xx[0]*ll;
+        double y = xx[1]*ll;
+        double z = (xx.Size()==3) ? xx[2]*ll : 0.0;
+
+        double r=std::sin(x)*std::cos(y) +
+                std::sin(y)*std::cos(z) +
+               std::sin(z)*std::cos(x) ;
+
+        if(r>0.0){return 1.0;}
+        return -1.0;
+    }
+
+private:
+    double ll;
+};
+
+
+class CheseCoeff:public Coefficient
+{
+public:
+    CheseCoeff(double cell_size=1.0)
+    {
+        ll=cell_size;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        double x = xx[0]*ll;
+        double y = xx[1]*ll;
+        double z = (xx.Size()==3) ? xx[2]*ll : 0.0;
+
+        double r=std::cos(x)*std::cos(y)*std::cos(z)-0.1;
+
+        if(r>0.0){return -1.0;}
+        return 1.0;
+
+    }
+
+private:
+    double ll;
+
+};
+
+class DispSol2D:public VectorCoefficient
+{
+public:
+    DispSol2D(double E_=1.0, double nu_=0.3,
+              double lx_=1.0,double ly_=1.0):VectorCoefficient(2)
+    {
+        E=E_;
+        nu=nu_;
+        lx=lx_;
+        ly=ly_;
+
+    }
+
+    virtual
+    void  Eval(Vector &V, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        V.SetSize(2);
+
+        V[0] = sin(lx*xx[0])+cos(ly*xx[1]);
+        V[1] = cos(lx*xx[0]+ly*xx[1]);
+    }
+
+private:
+
+    double E;
+    double nu;
+    double lx;
+    double ly;
+};
+
+class StressSol2D:public MatrixCoefficient
+{
+public:
+    StressSol2D(double E_=1.0, double nu_=0.3,
+                double lx_=1.0,double ly_=1.0):MatrixCoefficient(2)
+    {
+        E=E_;
+        nu=nu_;
+        lx=lx_;
+        ly=ly_;
+    }
+
+    virtual
+    void  Eval(DenseMatrix &ss, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+
+        ss(0,0) = 2.0*E/(2.0+2.0*nu)*lx*cos(lx*xx[0])
+                +E*nu/(1.0+nu)/(1.0-2.0*nu)*(lx*cos(lx*xx[0])-ly*sin(lx*xx[0]+ly*xx[1]));
+        ss(0,1) = -E*(ly*sin(ly*xx[1])+lx*sin(lx*xx[0]+ly*xx[1]))/(2.0+2.0*nu);
+        ss(1,0) = -E*(ly*sin(ly*xx[1])+lx*sin(lx*xx[0]+ly*xx[1]))/(2.0+2.0*nu);
+        ss(1,1) = -2.0*E/(2.0+2.0*nu)*ly*sin(lx*xx[0]+ly*xx[1])
+                +E*nu/(1.0+nu)/(1.0-2.0*nu)*(lx*cos(lx*xx[0])-ly*sin(lx*xx[0]+ly*xx[1]));
+    }
+
+private:
+
+    double E;
+    double nu;
+    double lx;
+    double ly;
+
+};
+
+class BdrLoadSol2D:public VectorCoefficient
+{
+public:
+    BdrLoadSol2D(StressSol2D* sco_, ParGridFunction* lsf_):VectorCoefficient(2)
+    {
+        sco=sco_;
+        lsf=lsf_;
+    }
+    virtual
+    void  Eval(Vector &vv, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        vv.SetSize(2);
+        Vector n(2);
+        DenseMatrix ss(2);
+        sco->Eval(ss,T,ip);
+        T.SetIntPoint(&ip);
+        lsf->GetGradient(T,n);
+        double nr=n.Norml2();
+        n/=-nr;
+        ss.Mult(n,vv);
+    }
+
+
+private:
+    StressSol2D* sco;
+    ParGridFunction* lsf;
+
+
+};
+
+class ForceSol2D:public VectorCoefficient
+{
+public:
+    ForceSol2D(double E_=1.0, double nu_=0.3,
+               double lx_=1.0,double ly_=1.0):VectorCoefficient(2)
+    {
+         E=E_;
+         nu=nu_;
+         lx=lx_;
+         ly=ly_;
+    }
+
+    virtual
+    void  Eval(Vector &ff, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        ff.SetSize(2);
+
+        ff[0] = 2.0*E/(2.0+2.0*nu)*lx*lx*sin(lx*xx[0])
+                -E*nu/(1.0+nu)/(1.0-2.0*nu)*(-lx*lx*sin(lx*xx[0])-ly*lx*cos(lx*xx[0]+ly*xx[1]))
+                +E*(ly*ly*cos(ly*xx[1])+ly*lx*cos(lx*xx[0]+ly*xx[1]))/(2.0+2.0*nu);
+
+        ff[1] = E*lx*lx*cos(lx*xx[0]+ly*xx[1])/(2.0+2.0*nu)
+                +2.0*E/(2.0+2.0*nu)*ly*ly*cos(lx*xx[0]+ly*xx[1])
+                +E*nu/(1.0+nu)/(1.0-2.0*nu)*ly*ly*cos(lx*xx[0]+ly*xx[1]);
+    }
+
+private:
+
+    double E;
+    double nu;
+    double lx;
+    double ly;
+};
+
+class DispSol3D:public VectorCoefficient
+{
+public:
+    DispSol3D(double E_=1.0, double nu_=0.3,
+              double lx_=1.0,double ly_=1.0, double lz_=1.0):VectorCoefficient(3)
+    {
+        E=E_;
+        nu=nu_;
+        lx=lx_;
+        ly=ly_;
+        lz=lz_;
+
+    }
+
+    virtual
+    void  Eval(Vector &u, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        u.SetSize(3);
+
+        u[0] = sin(lx*xx[0])+cos(ly*xx[1]);
+        u[1] = sin(lx*xx[1])+cos(lz*xx[2]);
+        u[2] = sin(lz*xx[2])+cos(lx*xx[0]);
+    }
+
+private:
+    double E;
+    double nu;
+    double lx;
+    double ly;
+    double lz;
+
+};
+
+
+
+
+class StressSol3D:public MatrixCoefficient
+{
+public:
+    StressSol3D(double E_=1.0, double nu_=0.3,
+              double lx_=1.0,double ly_=1.0, double lz_=1.0):MatrixCoefficient(3)
+    {
+        E=E_;
+        nu=nu_;
+        lx=lx_;
+        ly=ly_;
+        lz=lz_;
+   }
+
+    virtual
+    void  Eval(DenseMatrix &ss, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        double nnu=nu;
+
+        ss(0,0) = 2.0*E/(2.0+2.0*nnu)*lx*cos(lx*xx[0])+E*nnu/(1.0+nnu)/(1.0-2.0*
+        nnu)*(lx*cos(lx*xx[0])+lx*cos(lx*xx[1])+lz*cos(lz*xx[2]));
+
+        ss(0,1) = -E/(2.0+2.0*nnu)*ly*sin(ly*xx[1]);
+        ss(0,2) = -E/(2.0+2.0*nnu)*lx*sin(lx*xx[0]);
+        ss(1,0) = -E/(2.0+2.0*nnu)*ly*sin(ly*xx[1]);
+        ss(1,1) = 2.0*E/(2.0+2.0*nnu)*lx*cos(lx*xx[1])+E*nnu/(1.0+nnu)/(1.0-2.0*
+        nnu)*(lx*cos(lx*xx[0])+lx*cos(lx*xx[1])+lz*cos(lz*xx[2]));
+
+        ss(1,2) = -E/(2.0+2.0*nnu)*lz*sin(lz*xx[2]);
+        ss(2,0) = -E/(2.0+2.0*nnu)*lx*sin(lx*xx[0]);
+        ss(2,1) = -E/(2.0+2.0*nnu)*lz*sin(lz*xx[2]);
+        ss(2,2) = 2.0*E/(2.0+2.0*nnu)*lz*cos(lz*xx[2])+E*nnu/(1.0+nnu)/(1.0-2.0*
+        nnu)*(lx*cos(lx*xx[0])+lx*cos(lx*xx[1])+lz*cos(lz*xx[2]));
+    }
+
+
+
+private:
+    double E;
+    double nu;
+    double lx;
+    double ly;
+    double lz;
+
+};
+
+class ForceSol3D:public VectorCoefficient
+{
+public:
+    ForceSol3D(double E_=1.0, double nu_=0.3,
+              double lx_=1.0,double ly_=1.0, double lz_=1.0):VectorCoefficient(3)
+    {
+        E=E_;
+        nu=nu_;
+        lx=lx_;
+        ly=ly_;
+        lz=lz_;
+
+    }
+
+    virtual
+    void  Eval(Vector &ff, ElementTransformation &T,
+               const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        ff.SetSize(3);
+
+        double nnu=nu;
+        ff[0] = 2.0*E/(2.0+2.0*nnu)*lx*lx*sin(lx*xx[0])+E*nnu/(1.0+nnu)/(1.0-2.0*
+  nnu)*lx*lx*sin(lx*xx[0])+E/(2.0+2.0*nnu)*ly*ly*cos(ly*xx[1]);
+        ff[1] = 2.0*E/(2.0+2.0*nnu)*lx*lx*sin(lx*xx[1])+E*nnu/(1.0+nnu)/(1.0-2.0*
+  nnu)*lx*lx*sin(lx*xx[1])+E/(2.0+2.0*nnu)*lz*lz*cos(lz*xx[2]);
+        ff[2] = E/(2.0+2.0*nnu)*lx*lx*cos(lx*xx[0])+2.0*E/(2.0+2.0*nnu)*lz*lz*sin
+  (lz*xx[2])+E*nnu/(1.0+nnu)/(1.0-2.0*nnu)*lz*lz*sin(lz*xx[2]);
+
+    }
+
+private:
+    double E;
+    double nu;
+    double lx;
+    double ly;
+    double lz;
+
+};
+
+class StressCompCoef:public Coefficient
+{
+public:
+    StressCompCoef(StressSol2D* sco_,LinIsoElasticityCoefficient* lco_)
+    {
+        sco2d=sco_;
+        lco=lco_;
+    }
+
+    StressCompCoef(StressSol3D* sco_,LinIsoElasticityCoefficient* lco_)
+    {
+        sco2d=nullptr;
+        sco3d=sco_;
+        lco=lco_;
+    }
+
+
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        DenseMatrix ss(3);
+        lco->EvalStress(ss, T, ip);
+
+        DenseMatrix ssc(3);
+        double r=0.0;
+        if(sco2d!=nullptr){
+            ssc=0.0;
+            sco2d->Eval(ssc,T,ip);
+            for(int i=0;i<2;i++){
+            for(int j=0;j<2;j++){
+                r=r+(ssc(i,j)-ss(i,j))*(ssc(i,j)-ss(i,j));
+            }}
+        }else{
+            sco3d->Eval(ssc,T,ip);
+            for(int i=0;i<3;i++){
+            for(int j=0;j<3;j++){
+                r=r+(ssc(i,j)-ss(i,j))*(ssc(i,j)-ss(i,j));
+            }}
+        }
+
+        return sqrt(r);
+
+
+    }
+
+private:
+    StressSol3D* sco3d;
+    StressSol2D* sco2d;
+    LinIsoElasticityCoefficient* lco;
+};
+
+
+void CompDisplH1Norm(VectorCoefficient& ut, ParGridFunction* uu,
+                       Array<int> marks, CutIntegrationRules& cut_int, double& L2err, double& H1err)
+{
+    ParMesh* mesh=uu->ParFESpace()->GetParMesh();
+    int order = uu->ParFESpace()->GetMaxElementOrder();
+    int dim= mesh->SpaceDimension();
+    H1_FECollection fec(order+3,dim);
+    ParFiniteElementSpace fes(mesh,&fec,dim);
+
+    ParGridFunction up(&fes);
+    up.ProjectCoefficient(ut);
+
+    //add the cuts
+    double cerr=0.0;
+    double derr=0.0;
+    double rerr=0.0;
+
+    {
+        ElementTransformation *trans;
+        Vector duu(dim);
+        Vector dtt(dim);
+        DenseMatrix dduu(dim);
+        DenseMatrix ddtt(dim);
+        double w;
+        const IntegrationRule* ir;
+        for(int i=0;i<fes.GetNE();i++)
+        {
+            const FiniteElement* el=fes.GetFE(i);
+            //get the element transformation
+            trans = fes.GetElementTransformation(i);
+
+            if(marks[i]==ElementMarker::SBElementType::INSIDE){
+                ir=&IntRules.Get(el->GetGeomType(), order+3);
+            }else
+            if(marks[i]==ElementMarker::SBElementType::CUT){
+                ir=cut_int.GetSurfIntegrationRule(i);
+            }else{
+                continue;
+            }
+
+            for(int j=0; j<ir->GetNPoints();j++){
+                const IntegrationPoint &ip = ir->IntPoint(j);
+                trans->SetIntPoint(&ip);
+                w=ip.weight * trans->Weight();
+
+                uu->GetVectorValue(*trans,ip,duu);
+                up.GetVectorValue(*trans,ip,dtt);
+
+                duu.Add(-1.0,dtt);
+                cerr=cerr+w*(duu*duu);
+
+                uu->GetVectorGradient(*trans,dduu);
+                up.GetVectorGradient(*trans,ddtt);
+
+                dduu.Add(-1.0,ddtt);
+                for(int pp=0;pp<dim;pp++){
+                for(int kk=0;kk<dim;kk++){
+                    derr=derr+w*dduu(kk,pp)*dduu(kk,pp);
+                }}
+            }
+        }
+
+    }
+
+    MPI_Reduce(&cerr,&rerr,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
+    cerr=sqrt(rerr); L2err=cerr;
+    rerr=0.0;
+    MPI_Reduce(&derr,&rerr,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
+    cerr=cerr+sqrt(rerr); H1err=cerr;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int rs_levels = 2;
+   int order = 2;
+   int cut_int_order = order;
+   const char *device_config = "cpu";
+   double stiff_ratio=1e-6;
+   bool visualization = true;
+
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&cut_int_order, "-co", "--corder",
+                  "Cut integration order");
+   args.AddOption(&stiff_ratio,"-sr", "--stiff_ratio",
+                  "Stiffness ratio");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0) { args.PrintOptions(cout); }
+
+   // Enable hardware devices such as GPUs, and programming models such as CUDA,
+   // OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myrank == 0) { device.Print(); }
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+   if(myrank==0){
+       std::cout<<"Num elements="<<mesh.GetNE()<<std::endl;
+   }
+
+
+
+   // MPI distribution.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   FilterSolver* filter=new FilterSolver(0.1,&pmesh,2);
+
+   ParFiniteElementSpace* dfes=filter->GetDesignFES();
+   ParFiniteElementSpace* ffes=filter->GetFilterFES();
+
+   ParGridFunction desgf(dfes); desgf=0.0;
+   ParGridFunction filgf(ffes); filgf=0.0;
+
+   {// project the coefficient and filter
+       //GyroidCoeff gc(6.0*M_PI);
+       CheseCoeff  gc(2.0*M_PI);
+       desgf.ProjectCoefficient(gc);
+       Vector tdes(dfes->GetTrueVSize()); tdes=0.0;
+       desgf.GetTrueDofs(tdes);
+       Vector tfil(ffes->GetTrueVSize()); tfil=0.0;
+       filter->Mult(tdes,tfil);
+       filgf.SetFromTrueDofs(tfil);
+   }
+
+   ElementMarker* elmark=new ElementMarker(pmesh,false,true);
+   elmark->SetLevelSetFunction(filgf);
+
+   Array<int> marks;
+   elmark->MarkElements(marks);
+   Array<int> ghost_penalty_marks;
+   elmark->MarkGhostPenaltyFaces(ghost_penalty_marks);
+
+
+   //define the cut integration rules
+   CutIntegrationRules cut_int(2*cut_int_order, filgf, marks);
+
+   for(int i=0;i<pmesh.GetNE();i++){
+       pmesh.SetAttribute(i,marks[i]);
+   }
+
+   DispSol2D dsol2d(1.0,0.3,2.0*M_PI,2.0*M_PI);
+   StressSol2D ssol2d(1.0,0.3,2.0*M_PI,2.0*M_PI);
+   DispSol3D dsol3d(1.0,0.3,2.0*M_PI,2.0*M_PI,2.0*M_PI);
+   StressSol3D ssol3D(1.0,0.3,2.0*M_PI,2.0*M_PI,2.0*M_PI);
+
+
+
+
+   CFElasticitySolver* elsolv=new CFElasticitySolver(&pmesh,order);
+   elsolv->SetGhostPenalty(0.001,ghost_penalty_marks);
+
+   Vector vf(dim); vf=0.0; vf(1)=0.0;
+   //VectorConstantCoefficient* ff=new VectorConstantCoefficient(vf);
+   LinIsoElasticityCoefficient* lec=new LinIsoElasticityCoefficient(1.0,0.3);
+   elsolv->SetLinearSolver(1e-12,1e-12,400);
+   elsolv->SetNewtonSolver(1e-10,1e-12,20,1);
+   if(dim==2){
+       ForceSol2D* fsol2d= new ForceSol2D(1.0,0.3,2.0*M_PI,2.0*M_PI);
+       BdrLoadSol2D* surf_load2d=new BdrLoadSol2D(&ssol2d,&filgf);
+       elsolv->AddMaterial(lec,fsol2d,surf_load2d);
+       //elsolv->AddMaterial(lec,fsol2d,nullptr);
+       //elsolv->AddMaterial(lec,nullptr,nullptr);
+   }else{//3D
+       ForceSol3D* fsol3d=new ForceSol3D(1.0,0.3,2.0*M_PI,2.0*M_PI,2.0*M_PI);
+       elsolv->AddMaterial(lec,fsol3d,nullptr);
+   }
+   //elsolv->AddDispBC(2,4,0.0);
+   if(dim==2){
+       elsolv->AddDispBC(1,dsol2d);
+       elsolv->AddDispBC(2,dsol2d);
+       elsolv->AddDispBC(3,dsol2d);
+       elsolv->AddDispBC(4,dsol2d);
+   }else{ //3D
+       elsolv->AddDispBC(2,dsol3d);
+   }
+   elsolv->SetLSF(filgf,marks, cut_int);
+   elsolv->SetStiffnessRatio(stiff_ratio);
+   elsolv->FSolve();
+   ParGridFunction& u=elsolv->GetDisplacements();
+
+   //chack the displacement and the stress coefficients
+   //displacement field
+   ParGridFunction ug(u);
+   L2_FECollection l2fec(4,dim);
+   ParFiniteElementSpace l2fes(&pmesh,&l2fec,1);
+   ParGridFunction errgf(&l2fes);
+   if(dim==2){
+       ug.ProjectCoefficient(dsol2d);
+       lec->SetDisplacementField(ug);
+       StressSol2D sco(1.0,0.3,2.0*M_PI,2.0*M_PI);
+       StressCompCoef scco(&sco,lec);
+       errgf.ProjectCoefficient(scco);
+
+   }else{
+       ug.ProjectCoefficient(dsol3d);
+       lec->SetDisplacementField(ug);
+       StressSol3D sco(1.0,0.3,2.0*M_PI,2.0*M_PI,2.0*M_PI);
+       StressCompCoef scco(&sco,lec);
+       errgf.ProjectCoefficient(scco);
+   }
+
+
+   double L2err;
+   double H1err;
+   if(dim==2){
+       CompDisplH1Norm(dsol2d,&u,marks,cut_int,L2err,H1err);
+   }else{
+       CompDisplH1Norm(dsol3d,&u,marks,cut_int,L2err,H1err);
+   }
+
+
+   if(myrank==0){
+       std::cout.setf( std::ios_base::scientific, std::ios_base::floatfield );
+       std::cout<<"L2 err = "<<L2err<<std::endl;
+       std::cout<<"H1 err = "<<H1err<<std::endl;
+   }
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.SetHighOrderOutput(true);
+   dacol.RegisterField("design", &desgf);
+   dacol.RegisterField("flter", &filgf);
+   dacol.RegisterField("disp",&u);
+   dacol.RegisterField("disp_sol",&ug);
+   dacol.RegisterField("err",&errgf);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+
+   delete elsolv;
+   delete elmark;
+   delete filter;
+
+   return 0;
+}

miniapps/shifted/cut_obj_test.cpp

@@ -0,0 +1,439 @@
+#include <fstream>
+#include <iostream>
+
+#include "mfem.hpp"
+#include "marking.hpp"
+#include "mtop_solvers.hpp"
+#include "mtop_filters.hpp"
+#include "shape_grad.hpp"
+
+
+using namespace mfem;
+using namespace std;
+
+class GyroidCoeff:public Coefficient
+{
+public:
+    GyroidCoeff(double cell_size=1.0){
+        ll=cell_size;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+        double x = xx[0]*ll;
+        double y = xx[1]*ll;
+        double z = (xx.Size()==3) ? xx[2]*ll : 0.0;
+
+        double r=std::sin(x)*std::cos(y) +
+                std::sin(y)*std::cos(z) +
+               std::sin(z)*std::cos(x) ;
+
+        return r;
+
+        //if(r>0.0){return 1.0;}
+        //return -1.0;
+    }
+
+private:
+    double ll;
+};
+
+
+
+class CheseCoeff:public Coefficient
+{
+public:
+    CheseCoeff(double cell_size=1.0)
+    {
+        ll=cell_size;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+
+        double x = xx[0]*ll;
+        double y = xx[1]*ll;
+        double z = (xx.Size()==3) ? xx[2]*ll : 0.0;
+
+        double r=std::cos(x)*std::cos(y)*std::cos(z)-0.1;
+
+        if(r>0.0){return 0.0;}
+        return 1.0;
+    }
+
+private:
+    double ll;
+
+};
+
+class ThresholdCoeff:public Coefficient
+{
+public:
+    ThresholdCoeff(Coefficient& coeff_,Coefficient& threshold_)
+    {
+        coeff=&coeff_;
+        threshold=&threshold_;
+        loc_thres=nullptr;
+        loc_coeff=nullptr;
+    }
+
+
+    ThresholdCoeff(Coefficient& coeff_,double thr)
+    {
+        coeff=&coeff_;
+        loc_thres=new ConstantCoefficient(thr);
+        threshold=loc_thres;
+        loc_coeff=nullptr;
+
+    }
+
+    ThresholdCoeff(GridFunction& coeff_,double thr)
+    {
+        loc_coeff=new GridFunctionCoefficient(&coeff_);
+        coeff=loc_coeff;
+        loc_thres=new ConstantCoefficient(thr);
+        threshold=loc_thres;
+
+    }
+
+    ~ThresholdCoeff()
+    {
+        delete loc_thres;
+        delete loc_coeff;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        double cc=coeff->Eval(T,ip);
+        double hh=threshold->Eval(T,ip);
+        if(cc>hh){return 1.0;}
+        else{return 0.0;}
+    }
+
+private:
+    Coefficient* coeff;
+    Coefficient* threshold;
+    Coefficient* loc_thres;
+    Coefficient* loc_coeff;
+};
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int solver_type = 0;
+   int rs_levels = 2;
+   int order = 2;
+   int cut_int_order = order;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   double stiff_ratio=1e-6;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&stiff_ratio,"-sr", "--stiff_ratio",
+                  "Stiffness ratio");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0) { args.PrintOptions(cout); }
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.SpaceDimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+   if(myrank==0){
+       std::cout<<"Num elements="<<mesh.GetNE()<<std::endl;
+   }
+
+   // MPI distribution.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   FilterSolver* filter=new FilterSolver(0.1,&pmesh,2);
+   filter->SetSolver(1e-8,1e-12,100,0);
+if(myrank==0){std::cout<<"step 1"<<std::endl;}
+
+   ParGridFunction desgf(filter->GetDesignFES());
+   ParGridFunction filgf(filter->GetFilterFES());
+
+   {
+       //project the coefficient and filter
+       GyroidCoeff gc(6.0*M_PI);
+       ThresholdCoeff tc(gc,0.0);
+
+       CheseCoeff cc(2*M_PI);
+       desgf.ProjectCoefficient(cc);
+       Vector tdes(filter->GetDesignFES()->GetTrueVSize()); tdes=0.0;
+       desgf.GetTrueDofs(tdes);
+       Vector tfil(filter->GetFilterFES()->GetTrueVSize()); tfil=0.0;
+       filter->Mult(tdes,tfil);
+       filgf.SetFromTrueDofs(tfil);
+   }
+
+if(myrank==0){std::cout<<"step 1"<<std::endl;}
+
+   ElementMarker* elmark=new ElementMarker(pmesh,false,true);
+   ParGridFunction lsgf(filgf); lsgf-=0.5;
+   elmark->SetLevelSetFunction(lsgf);
+   Array<int> marks;
+   elmark->MarkElements(marks);
+
+if(myrank==0){std::cout<<"step 2"<<std::endl;}
+   //define the cut integration rule
+   CutIntegrationRules* cut_int=new CutIntegrationRules(2*cut_int_order, lsgf,marks);
+if(myrank==0){std::cout<<"step 3"<<std::endl;}
+   CFElasticitySolver* elsolv=new CFElasticitySolver(&pmesh,order);
+   Vector vf(dim); vf=0.0; vf(1)=1.0;
+   VectorConstantCoefficient* ff=new VectorConstantCoefficient(vf);
+   LinIsoElasticityCoefficient* lec=new LinIsoElasticityCoefficient(1.0,0.3);
+   elsolv->AddMaterial(lec,ff,nullptr);
+   elsolv->SetLinearSolver(1e-12,1e-12,400);
+   elsolv->SetNewtonSolver(1e-10,1e-12,20,0);
+   elsolv->SetLSF(lsgf,marks,*cut_int);
+   elsolv->AddDispBC(2,4,0.0);
+   elsolv->SetStiffnessRatio(stiff_ratio);
+   elsolv->FSolve();
+   ParGridFunction& u=elsolv->GetDisplacements();
+
+
+if(myrank==0){std::cout<<"step 4"<<std::endl;}
+   ComplianceObjectiveCut* cobj=new ComplianceObjectiveCut();
+   cobj->SetCutIntegrationRules(marks,cut_int);
+   cobj->SetE(1.0);
+   cobj->SetPoissonRatio(0.3);
+   cobj->SetVolForce(ff);
+   double cpl=cobj->Eval(u,lsgf);
+   std::cout<<"Compliance="<<cpl<<std::endl;
+   Vector cfgrad(filter->GetFilterFES()->GetTrueVSize()); cfgrad=0.0;
+   Vector cograd(filter->GetDesignFES()->GetTrueVSize()); cograd=0.0;
+
+   cobj->Grad(u,lsgf,cfgrad);
+   filter->MultTranspose(cfgrad,cograd);
+
+   ParGridFunction cgr(filter->GetFilterFES()); cgr.SetFromTrueDofs(cfgrad);
+   ParGridFunction cor(filter->GetDesignFES()); cor.SetFromTrueDofs(cograd);
+
+
+   //displ objective
+   DisplObjectiveCut* dobj=new DisplObjectiveCut();
+   GyroidCoeff gc(2.0);
+   GyroidCoeff gu(4.0);
+   VectorArrayCoefficient utg(dim);
+   for(int i=0;i<dim;i++){
+       utg.Set(i,&gu,false);
+   }
+   dobj->SetCutIntegrationRules(cut_int);
+   dobj->SetE(1.0);
+   dobj->SetPoissonRatio(0.3);
+   dobj->SetVolForce(ff);
+   dobj->SetPower(2.0);
+   dobj->SetSolver(elsolv);
+   dobj->SetIntegrationWeight(&gc);
+   dobj->SetTargetDisplacement(&utg);
+
+   double dpl=dobj->Eval(lsgf);
+   if(myrank==0){
+       std::cout<<" Displ obj="<<dpl<<std::endl;
+   }
+
+   Vector dfgrad(filter->GetFilterFES()->GetTrueVSize()); dfgrad=0.0;
+   Vector dograd(filter->GetDesignFES()->GetTrueVSize()); dograd=0.0;
+
+   std::cout<<"grad in"<<std::endl;
+   dobj->Grad(lsgf,dfgrad);
+   std::cout<<"grad on"<<std::endl;
+   filter->MultTranspose(dfgrad,dograd);
+   std::cout<<"grad out"<<std::endl;
+
+   ParGridFunction dgr(filter->GetFilterFES()); dgr.SetFromTrueDofs(dfgrad);
+   ParGridFunction dor(filter->GetDesignFES()); dor.SetFromTrueDofs(dograd);
+
+
+   //stress objective
+   StressObjectiveCut* sobj=new StressObjectiveCut();
+   sobj->SetCutIntegrationRules(cut_int);
+   sobj->SetE(1.0);
+   sobj->SetPoissonRatio(0.3);
+   sobj->SetVolForce(ff);
+   sobj->SetPower(2.0);
+   sobj->SetSolver(elsolv);
+   sobj->SetIntegrationWeight(&gc);
+
+   double spl=sobj->Eval(lsgf);
+   if(myrank==0){
+       std::cout<<" Stress obj="<<spl<<std::endl;
+   }
+
+   Vector sfgrad(filter->GetFilterFES()->GetTrueVSize()); sfgrad=0.0;
+   Vector sograd(filter->GetDesignFES()->GetTrueVSize()); sograd=0.0;
+
+   std::cout<<"grad in"<<std::endl;
+   sobj->Grad(lsgf,sfgrad);
+   std::cout<<"grad on"<<std::endl;
+   filter->MultTranspose(sfgrad,sograd);
+   std::cout<<"grad out"<<std::endl;
+
+   ParGridFunction sgr(filter->GetFilterFES()); sgr.SetFromTrueDofs(sfgrad);
+   ParGridFunction sor(filter->GetDesignFES()); sor.SetFromTrueDofs(sograd);
+
+
+   VolObjectiveCut* vobj=new VolObjectiveCut();
+   vobj->SetCutIntegrationRules(marks,cut_int);
+
+
+   double vol00=vobj->Eval(lsgf);
+   Vector fvolgrad(filter->GetFilterFES()->GetTrueVSize());
+   Vector ovolgrad(filter->GetDesignFES()->GetTrueVSize());
+   vobj->Grad(lsgf,fvolgrad);
+   filter->MultTranspose(fvolgrad,ovolgrad);
+   ParGridFunction vgr(filter->GetFilterFES()); vgr.SetFromTrueDofs(fvolgrad);
+   ParGridFunction vor(filter->GetDesignFES()); vor.SetFromTrueDofs(ovolgrad);
+
+   std::cout<<"vol="<<vol00<<std::endl;
+
+   //FD check
+   {
+       mfem::Vector prtv;
+       mfem::Vector tmpv;
+
+       mfem::Vector lsfv; desgf.GetTrueDofs(lsfv);
+
+       prtv.SetSize(ovolgrad.Size());
+       tmpv.SetSize(ovolgrad.Size());
+
+       prtv.Randomize();
+       double nd=mfem::InnerProduct(pmesh.GetComm(),prtv,prtv);
+       double td=mfem::InnerProduct(pmesh.GetComm(),prtv,ovolgrad);
+       double cd=mfem::InnerProduct(pmesh.GetComm(),prtv,cograd);
+       double dd=mfem::InnerProduct(pmesh.GetComm(),prtv,dograd);
+       double ss=mfem::InnerProduct(pmesh.GetComm(),prtv,sograd);
+
+       td=td/nd;
+       cd=cd/nd;
+       dd=dd/nd;
+       ss=ss/nd;
+       double lsc=1.0;
+       double lqoi;
+       double dqoi;
+       double sqoi;
+
+       for(int l=0; l<8;l++){
+           lsc/=10.0;
+           prtv/=10.0;
+
+           add(prtv,lsfv,tmpv);
+           filter->Mult(tmpv,fvolgrad);
+           lsgf.SetFromTrueDofs(fvolgrad); lsgf-=0.5;
+           elmark->SetLevelSetFunction(lsgf);
+           elmark->MarkElements(marks);
+           delete cut_int;
+           cut_int=new CutIntegrationRules(2*cut_int_order, lsgf, marks);
+
+           //compute the linear elastic model
+           elsolv->SetLSF(lsgf,marks,*cut_int);
+           elsolv->FSolve();
+           u=elsolv->GetDisplacements();
+
+           cobj->SetCutIntegrationRules(marks,cut_int);
+           double cqoi=cobj->Eval(u,lsgf);
+           double sd=(cqoi-cpl)/lsc;
+
+           dobj->SetCutIntegrationRules(cut_int);
+           dqoi=dobj->Eval(lsgf);
+           double md=(dqoi-dpl)/lsc;
+
+           sobj->SetCutIntegrationRules(cut_int);
+           sqoi=sobj->Eval(lsgf);
+           double ssd=(sqoi-spl)/lsc;
+
+
+
+           vobj->SetCutIntegrationRules(marks, cut_int);
+           lqoi=vobj->Eval(lsgf);
+           double ld=(lqoi-vol00)/lsc;
+           if(myrank==0){
+               std::cout<<" obj="<<vol00<<" lvo="<< lqoi<<" dx="<<lsc<<" FD app="<< ld/nd<<
+                          " gr="<< td <<" err="<< std::fabs(ld/nd-td) <<std::endl;
+
+               std::cout<<" cpl="<<cpl<<" lcpl="<< cqoi<<" dx="<<lsc<<" FD app="<< sd/nd<<
+                          " gr="<< cd <<" err="<< std::fabs(sd/nd-cd) <<std::endl;
+
+               std::cout<<" dpl="<<dpl<<" lcpl="<< dqoi<<" dx="<<lsc<<" FD app="<< md/nd<<
+                          " gr="<< dd <<" err="<< std::fabs(md/nd-dd) <<std::endl;
+
+               std::cout<<" spl="<<spl<<" scpl="<< sqoi<<" dx="<<lsc<<" FD app="<< ssd/nd<<
+                          " gr="<< ss <<" err="<< std::fabs(ssd/nd-ss) <<std::endl;
+           }
+       }
+   }
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.SetHighOrderOutput(true);
+   dacol.RegisterField("design", &desgf);
+   dacol.RegisterField("flter", &filgf);
+   dacol.RegisterField("fgrad",&vgr);
+   dacol.RegisterField("ograd",&vor);
+   dacol.RegisterField("lsf",&lsgf);
+   dacol.RegisterField("disp",&u);
+   dacol.RegisterField("cfgrad",&cgr);
+   dacol.RegisterField("cograd",&cor);
+   dacol.RegisterField("dfgrad",&dgr);
+   dacol.RegisterField("dograd",&dor);
+   dacol.RegisterField("adjoint",&(elsolv->GetADisplacements()));
+   dacol.RegisterField("sfgrad",&sgr);
+   dacol.RegisterField("sograd",&sor);
+
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+
+   delete cobj;
+   delete sobj;
+   delete dobj;
+   delete elsolv;
+   delete vobj;
+   delete cut_int;
+   delete elmark;
+   delete filter;
+}

miniapps/shifted/exo2gmsh.cpp

@@ -0,0 +1,151 @@
+#include "mfem.hpp"
+
+#ifdef MFEM_USE_PETSC
+#include "exodusII.h"
+#endif
+
+int main(int argc, char *argv[])
+{
+    // Parse command-line options.
+    const char *mesh_file = "input.exo";
+    const char *gmsh_file = "gmsh.out";
+
+
+    mfem::OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh",
+                   "Mesh file to use.");
+
+    args.AddOption(&gmsh_file, "-o", "--output",
+                   "GMSH output file.");
+
+
+    args.Parse();
+    if (!args.Good())
+    {
+       args.PrintUsage(std::cout);
+       return 1;
+    }
+    args.PrintOptions(std::cout);
+#ifdef MFEM_USE_PETSC
+
+    std::fstream gmsh;
+    gmsh.open(gmsh_file,std::ios::out);
+    gmsh<<"$MeshFormat\n2 0 8\n$EndMeshFormat"<<std::endl;
+
+    int comp_ws=sizeof(double);
+    int io_ws=0;
+    float version;
+
+    int exo_id=ex_open(mesh_file,EX_READ,&comp_ws,&io_ws,&version);
+
+    int num_dim, num_nodes, num_elems;
+    int num_elem_blk, num_node_sets, num_side_sets;
+    char title[256];
+
+
+    ex_get_init(exo_id, title, &num_dim, &num_nodes, &num_elems, &num_elem_blk, &num_node_sets,
+                      &num_side_sets);
+
+
+    std::cout<<"num_dim="<<num_dim<<std::endl;
+    std::cout<<"num_nodes="<<num_nodes<<std::endl;
+    std::cout<<"num_elems="<<num_elems<<std::endl;
+    std::cout<<"num_elem_blk="<<num_elem_blk<<std::endl;
+    std::cout<<"num_node_sets="<<num_node_sets<<std::endl;
+    std::cout<<"num_side_sets="<<num_side_sets<<std::endl;
+
+    //read the nodes and dump them to a file
+    {
+        int numn=num_nodes;
+        //ex_inquire(exo_id,EX_INQ_NODES,&numn ,NULL,NULL);
+        std::cout<<"Num nodes ="<<numn<<std::endl;
+
+        double* x=new double[numn];
+        double* y=new double[numn];
+        double* z=new double[numn];
+
+        ex_get_coord(exo_id,x,y,z);
+        gmsh<<"$Nodes"<<std::endl;
+        gmsh<<numn<<std::endl;
+
+        for(int i=0;i<numn;i++)
+        {
+            gmsh<<i+1<<" "<<x[i]<<" "<<y[i]<<" "<<z[i]<<std::endl;
+        }
+
+        gmsh<<"$EndNodes"<<std::endl;
+
+        delete [] x;
+        delete [] y;
+        delete [] z;
+    }
+
+    
+    //do the element processing
+    {
+        std::cout<<"Num blocks="<<num_elem_blk<<std::endl;
+        int ids[num_elem_blk];
+        ex_get_ids(exo_id, EX_ELEM_BLOCK, ids);
+        char el_type[MAX_STR_LENGTH+1];
+        int numel;
+        int numel_nod;
+        int numel_attr;
+
+        gmsh<<"$Elements"<<std::endl;
+        gmsh<<num_elems<<std::endl;
+
+        int el_id=0;
+
+        for(int bl=0;bl<num_elem_blk;bl++)
+        {
+            ex_get_block(exo_id,EX_ELEM_BLOCK,ids[bl],
+                         el_type,&numel,&numel_nod,nullptr,nullptr,&numel_attr);
+
+            std::cout<<"bl="<<ids[bl]<<" nel="<<numel;
+            std::cout<<" nod="<<numel_nod<<" natt="<<numel_attr;
+            std::cout<<std::endl;
+            int* conn=new int[numel*numel_nod];
+
+            ex_get_conn(exo_id,EX_ELEM_BLOCK,ids[bl],conn,nullptr,nullptr);
+
+            for(int i=0;i<numel;i++)
+            {
+                el_id++;
+                gmsh<<el_id<<" ";
+                //element type
+                if(strcasecmp(el_type,"TETRA4")==0){
+                    gmsh<<"4 ";
+                }else if(strcasecmp(el_type,"TETRA10")==0){
+                    gmsh<<"11 ";
+                }
+                else if(strcasecmp(el_type,"HEX8")==0){
+                    gmsh<<"5 ";
+                }else if(strcasecmp(el_type,"BAR2")==0){
+                    gmsh<<"1 ";
+                }else if(strcasecmp(el_type,"TRI3")==0){
+                    gmsh<<"2 ";
+                }else if(strcasecmp(el_type,"SHELL4")==0){
+                    gmsh<<"3 ";
+                }
+
+                //dump the rest of the element data
+                gmsh<<"3 "<<ids[bl]<<" "<<ids[bl];
+                gmsh<<" 0"; //num tags, geom domain, volume
+
+                for(int ii=0;ii<numel_nod;ii++)
+                {
+                    gmsh<<" "<<conn[i*numel_nod+ii];
+                }
+                gmsh<<std::endl;
+              }
+
+              delete [] conn;
+        }
+        gmsh<<"$EndElements"<<std::endl;
+    }
+    
+
+    gmsh.close();
+    ex_close(exo_id);
+#endif
+}

miniapps/shifted/marking.cpp

@@ -11,9 +11,388 @@
 
 #include "marking.hpp"
 
+#include "integ_algoim.hpp"
+
 namespace mfem
 {
 
+void ElementMarker::SetLevelSetFunction(Coefficient &ls_fun)
+{
+    FiniteElementCollection* fec=new H1_FECollection(h1_order,pmesh->Dimension());
+    ParFiniteElementSpace* pfes_sltn=new ParFiniteElementSpace(pmesh,fec);
+
+    Vector vals;
+    Array<int> vdofs;
+
+    if(use_cut_marks==false){
+    if(include_cut_elements){
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            {
+                int n = ir.GetNPoints();
+                vals.SetSize(n);
+                ElementTransformation *Tr = pfes_sltn->GetElementTransformation(e);
+                for(int k=0;k<n;k++){
+                    Tr->SetIntPoint(&ir.IntPoint(k));
+                    vals[k]=ls_fun.Eval(*Tr,ir.IntPoint(k));
+                }
+            }
+
+
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) { countp++; }
+                else{countn++;}
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }
+        }
+    }else{//DEFAULT - do not include cuts
+        elgf=(double)(SBElementType::OUTSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            {
+                int n = ir.GetNPoints();
+                vals.SetSize(n);
+                ElementTransformation *Tr = pfes_sltn->GetElementTransformation(e);
+                for(int k=0;k<n;k++){
+                    Tr->SetIntPoint(&ir.IntPoint(k));
+                    vals[k]=ls_fun.Eval(*Tr,ir.IntPoint(k));
+                }
+            }
+
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) { countp++; }
+                else{countn++;}
+            }
+            if (countp == ir.GetNPoints()) // completely inside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::INSIDE;
+            }
+        }
+    }}else{// use CUT mark
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            {
+                int n = ir.GetNPoints();
+                vals.SetSize(n);
+                ElementTransformation *Tr = pfes_sltn->GetElementTransformation(e);
+                for(int k=0;k<n;k++){
+                    Tr->SetIntPoint(&ir.IntPoint(k));
+                    vals[k]=ls_fun.Eval(*Tr,ir.IntPoint(k));
+                }
+            }
+
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0) {countp++;}
+                else {countn++;}
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }else
+            if ((countp>0)&&(countn>0))
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::CUT;
+            }
+        }
+    }
+    elgf.ExchangeFaceNbrData();
+
+    delete pfes_sltn;
+    delete fec;
+
+}
+
+void ElementMarker::SetLevelSetFunction(const ParGridFunction &ls_fun)
+{
+    ParFiniteElementSpace* pfes_sltn=ls_fun.ParFESpace();
+    Vector vals;
+    Array<int> vdofs;
+
+    if(use_cut_marks==false){
+    if(include_cut_elements){
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            ls_fun.GetValues(e, ir, vals);
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) {}
+                else{ countn++; }
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }
+        }
+    }else{//DEFAULT - do not include cuts
+        elgf=(double)(SBElementType::OUTSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            ls_fun.GetValues(e, ir, vals);
+            int countp = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) { countp++; }
+            }
+            if (countp == ir.GetNPoints()) // completely inside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::INSIDE;
+            }
+        }
+    }}else{//use CUT marks
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            ls_fun.GetValues(e, ir, vals);
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0) {countp++;}
+                else {countn++;}
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }else
+            if ((countp>0)&&(countn>0))
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::CUT;
+            }
+        }
+    }
+    if(pmesh->GetNRanks()>0){
+    elgf.ExchangeFaceNbrData();}
+
+}
+
+void ElementMarker::MarkElements(Array<int> &elem_marker)
+{
+    elem_marker.SetSize(pmesh->GetNE());
+    for(int e=0;e<pmesh->GetNE();e++)
+    {
+        ElementTransformation* tr=elfes->GetElementTransformation(e);
+        IntegrationPoint ip; ip.Init(0);
+        elem_marker[e] = elgf.GetValue(*tr, ip);
+    }
+}
+
+void ElementMarker::MarkGhostPenaltyFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(pmesh->GetNumFaces());
+    face_marker=SBFaceType::UNDEFINED;
+    IntegrationPoint ip; ip.Init(0);
+
+    for(int f=0;f<pmesh->GetNumFaces();f++){
+        auto *ft = pmesh->GetFaceElementTransformations(f, 3);
+        if (ft->Elem2No < 0) { continue; } //do not mark boundary faces
+        const int attr1 = elgf.GetValue(*ft->Elem1,ip);
+        const int attr2 = elgf.GetValue(*ft->Elem2,ip);
+        if((attr1==SBElementType::CUT)&&(attr2!=SBElementType::OUTSIDE))
+        {
+            face_marker[f]=SBFaceType::GHOSTP;
+        }else
+        if((attr1!=SBElementType::OUTSIDE)&&(attr2==SBElementType::CUT))
+        {
+            face_marker[f]=SBFaceType::GHOSTP;
+        }
+    }
+
+    elgf.ExchangeFaceNbrData();
+
+    for (int f = 0; f < pmesh->GetNSharedFaces(); f++)
+    {
+        auto *ftr = pmesh->GetSharedFaceTransformations(f, true);
+        const int attr1 = elgf.GetValue(*ftr->Elem1, ip);
+        const int attr2 = elgf.GetValue(*ftr->Elem2, ip);
+        int faceno = pmesh->GetSharedFace(f);
+        if((attr1==SBElementType::CUT)&&(attr2!=SBElementType::OUTSIDE))
+        {
+            face_marker[faceno]=SBFaceType::GHOSTP;
+        }else
+        if((attr1!=SBElementType::OUTSIDE)&&(attr2==SBElementType::CUT))
+        {
+            face_marker[faceno]=SBFaceType::GHOSTP;
+        }
+    }
+}
+
+
+void ElementMarker::MarkFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(pmesh->GetNumFaces());
+    face_marker=SBFaceType::UNDEFINED;
+    IntegrationPoint ip; ip.Init(0);
+
+    if(include_cut_elements==true){
+        for(int f=0;f<pmesh->GetNumFaces();f++){
+            auto *ft = pmesh->GetFaceElementTransformations(f, 3);
+            if (ft->Elem2No < 0) { continue; } //do not mark boundary faces
+            const int attr1 = elgf.GetValue(*ft->Elem1,ip);
+            const int attr2 = elgf.GetValue(*ft->Elem2,ip);
+            if((attr1==SBElementType::OUTSIDE)||(attr2==SBElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+    }else{
+        for(int f=0;f<pmesh->GetNumFaces();f++){
+            auto *ft = pmesh->GetFaceElementTransformations(f, 3);
+            if (ft->Elem2No < 0) { continue; } //do not mark boundary faces
+            const int attr1 = elgf.GetValue(*ft->Elem1,ip);
+            const int attr2 = elgf.GetValue(*ft->Elem2,ip);
+            if((attr1==SBElementType::INSIDE)||(attr2==SBElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+    }
+
+    elgf.ExchangeFaceNbrData();
+    if(include_cut_elements==true){
+        for (int f = 0; f < pmesh->GetNSharedFaces(); f++)
+        {
+            auto *ftr = pmesh->GetSharedFaceTransformations(f, true);
+            const int attr1 = elgf.GetValue(*ftr->Elem1, ip);
+            const int attr2 = elgf.GetValue(*ftr->Elem2, ip);
+            int faceno = pmesh->GetSharedFace(f);
+            if((attr1==SBElementType::OUTSIDE)||(attr2==SBElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[faceno]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+    }else{
+        for (int f = 0; f < pmesh->GetNSharedFaces(); f++)
+        {
+            auto *ftr = pmesh->GetSharedFaceTransformations(f, true);
+            const int attr1 = elgf.GetValue(*ftr->Elem1, ip);
+            const int attr2 = elgf.GetValue(*ftr->Elem2, ip);
+            int faceno = pmesh->GetSharedFace(f);
+            if((attr1==SBElementType::INSIDE)||(attr2==SBElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[faceno]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+
+    }
+}
+
+void ElementMarker::ListEssentialTDofs(const Array<int> &elem_marker,
+                                       ParFiniteElementSpace &lfes,
+                                       Array<int> &ess_tdof_list) const
+{
+    Array<int> dofs;
+
+    mfem::Vector vvdof; vvdof.SetSize(lfes.GetVSize()); vvdof=0.0;
+
+    for(int i=0;i<lfes.GetNE();i++)
+    {
+        if(elem_marker[i]==SBElementType::INSIDE){
+            lfes.GetElementVDofs(i,dofs);
+            for(int j=0;j<dofs.Size();j++){
+                vvdof[dofs[j]]=1.0;
+            }
+        }
+
+        if(include_cut_elements==true){
+            if(elem_marker[i]==SBElementType::CUT){
+                lfes.GetElementVDofs(i,dofs);
+                for(int j=0;j<dofs.Size();j++){
+                    vvdof[dofs[j]]=1.0;
+                }
+            }
+        }
+    }
+
+    Array<int> tdof_mark; tdof_mark.SetSize(lfes.GetTrueVSize());
+    Vector vtdof; vtdof.SetSize(lfes.GetTrueVSize()); vtdof=0.0;
+    lfes.GetProlongationMatrix()->MultTranspose(vvdof,vtdof);
+    for(int i=0;i<vtdof.Size();i++){
+        if(vtdof[i]<1.0){tdof_mark[i]=1;}
+        else{tdof_mark[i]=0;}
+    }
+
+    lfes.MarkerToList(tdof_mark, ess_tdof_list);
+}
+
+#ifdef MFEM_USE_ALGOIM
+
+
+CutIntegrationRules::CutIntegrationRules(int int_order,
+                                         ParGridFunction& lsf,
+                                         Array<int>& elm_markers):markers(elm_markers)
+{
+
+     AlgoimIntegrationRule* air;
+     Vector vlsf;//vector for the level-set-function
+     //cut element
+     Array<int> vdofs;
+     FiniteElementSpace* fespace=lsf.FESpace();
+     ElementTransformation *Tr;
+
+     vir.SetSize(fespace->GetNE()); vir=nullptr;
+     sir.SetSize(fespace->GetNE()); sir=nullptr;
+
+     for (int i=0; i<fespace->GetNE(); i++){
+     if(elm_markers[i]==ElementMarker::SBElementType::CUT){
+         const FiniteElement* le=fespace->GetFE(i);
+         fespace->GetElementDofs(i,vdofs);
+         Tr=fespace->GetElementTransformation(i);
+         lsf.GetSubVector(vdofs,vlsf);
+         air=new AlgoimIntegrationRule(int_order,*le,*Tr,vlsf);
+
+         vir[i]=new IntegrationRule(*(air->GetVolumeIntegrationRule()));
+         sir[i]=new IntegrationRule(*(air->GetSurfaceIntegrationRule()));
+
+         for (int j = 0; j < sir[i]->GetNPoints(); j++){
+              IntegrationPoint &ip = sir[i]->IntPoint(j);
+              Tr->SetIntPoint(&ip);
+              if(le->GetDim()==2){
+                  ip.weight=ip.weight*sqrt(Tr->Weight());
+              }else{
+                  ip.weight=ip.weight*pow(Tr->Weight(), 2.0/3.0);
+              }
+         }
+
+         delete air;
+     }}
+
+}
+
+CutIntegrationRules::~CutIntegrationRules()
+{
+    for(int i=0;i<vir.Size();i++)
+    {
+        delete vir[i];
+        delete sir[i];
+    }
+}
+
+
+#endif
+
+
 void ShiftedFaceMarker::MarkElements(const ParGridFunction &ls_func,
                                      Array<int> &elem_marker)
 {

miniapps/shifted/marking.hpp

@@ -16,6 +16,124 @@
 
 namespace mfem
 {
+// Marking operations for elements, faces, dofs, etc, related to shifted
+// boundary and interface methods.
+class ElementMarker{
+public:
+    enum SBElementType {INSIDE = 0, OUTSIDE = 1, CUT = 2};
+
+    enum SBFaceType {UNDEFINED = 0, SURROGATE = 1, GHOSTP = 2};
+
+    ///Defines element marker class with options to include the cut elements
+    /// (include_cut=true) or to mark the cut elements as SBElementType::CUT.
+    /// If use_cut=false the marking will use only INSIDE/OUTSIDE marks.
+    /// The last integer argument determines the order of the surrogate H1 field
+    /// for checking if an element is cut by a zero level set of an implicit
+    /// material distribution.
+    ElementMarker(ParMesh& mesh, bool include_cut=false,
+                  bool use_cut=false, int h1_order_=2)
+    {
+        pmesh=&mesh;
+        const int dim=pmesh->SpaceDimension();
+        elfec=new L2_FECollection(0,dim);
+        elfes=new ParFiniteElementSpace(pmesh,elfec,1);
+        elgf.SetSpace(elfes);
+        include_cut_elements=include_cut;
+        use_cut_marks=use_cut;
+
+        h1_order=h1_order_;
+    }
+
+    /// Destructor of the ElementMarker class
+    ~ElementMarker()
+    {
+        delete elfes;
+        delete elfec;
+    }
+
+    /// Mark elements according to the specified level-set
+    /// function.
+    void SetLevelSetFunction(const ParGridFunction& ls_fun);
+
+    /// Mark the elements according to the specified coefficient.
+    void SetLevelSetFunction(Coefficient& ls_fun);
+
+    /// Returns the marking of all the elements
+    /// in the mesh using the @a SBElementType
+    void MarkElements(Array<int> &elem_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a SBFaceType
+    void MarkFaces(Array<int> &face_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a SBFaceType.
+    /// The marks of all cut and faces between
+    /// cut and inside elements are set to GHOSTP
+    void MarkGhostPenaltyFaces(Array<int> &face_marker);
+
+    /// Lists all inactive dofs, i.e.,
+    ///  all dofs in the outside region.
+    void ListEssentialTDofs(const Array<int> &elem_marker,
+                            ParFiniteElementSpace &lfes,
+                            Array<int> &ess_tdof_list) const;
+private:
+    ParMesh* pmesh;
+    FiniteElementCollection* elfec;
+    ParFiniteElementSpace* elfes;
+    ParGridFunction elgf;
+
+    bool include_cut_elements;
+    bool use_cut_marks;
+
+    int h1_order; //order of the H1 FE space for level set functions defined by coefficient
+};
+
+
+#ifdef MFEM_USE_ALGOIM
+class CutIntegrationRules
+{
+public:
+    CutIntegrationRules(int int_order, ParGridFunction& lsf, Array<int>& elm_markers);
+    ~CutIntegrationRules();
+
+    ///Returns volumentric integration rulles for the cut elements.
+    ///The integration rule is different than a null pointer only
+    /// for cut elements.
+    const Array<IntegrationRule*>* GetVolIntegrationRule()
+    {
+        return &vir;
+    }
+
+    ///Returns the volumetric integration rule for element el.
+    const IntegrationRule* GetVolIntegrationRule(int el){
+        return vir[el];
+    }
+
+    ElementMarker::SBElementType GetElementMarker(int el){
+        return ElementMarker::SBElementType(markers[el]);
+    }
+
+    ///Returns surface integration rulles for the cut elements.
+    ///The integration rule is different than a null pointer only
+    /// for cut elements.
+    const Array<IntegrationRule*>* GetSurfIntegrationRule()
+    {
+        return &sir;
+    }
+
+    ///Returns surface integration rule for element el
+    const IntegrationRule* GetSurfIntegrationRule(int el)
+    {
+        return sir[el];
+    }
+
+private:
+    Array<IntegrationRule*> vir;
+    Array<IntegrationRule*> sir;
+    Array<int>& markers;
+};
+#endif
 
 // Marking operations for elements, faces, dofs, etc, related to shifted
 // boundary and interface methods.

miniapps/shifted/mtop_filters.cpp

@@ -0,0 +1 @@
+#include "mtop_filters.hpp"

miniapps/shifted/mtop_filters.hpp

@@ -0,0 +1,657 @@
+#ifndef MTOP_FILTERS_HPP
+#define MTOP_FILTERS_HPP
+
+#include "mfem.hpp"
+
+namespace mfem {
+
+namespace PointwiseTrans
+{
+
+/*  Standrd "Heaviside" projection in topology optimization with threshold eta
+ * and steepness of the projection beta.
+ * */
+inline
+double HProject(double rho, double eta, double beta)
+{
+    // tanh projection - Wang&Lazarov&Sigmund2011
+    double a=std::tanh(eta*beta);
+    double b=std::tanh(beta*(1.0-eta));
+    double c=std::tanh(beta*(rho-eta));
+    double rez=(a+c)/(a+b);
+    return rez;
+}
+
+/// Gradient of the "Heaviside" projection with respect to rho.
+inline
+double HGrad(double rho, double eta, double beta)
+
+{
+    double c=std::tanh(beta*(rho-eta));
+    double a=std::tanh(eta*beta);
+    double b=std::tanh(beta*(1.0-eta));
+    double rez=beta*(1.0-c*c)/(a+b);
+    return rez;
+}
+
+/// Second derivative of the "Heaviside" projection with respect to rho.
+inline
+double HHess(double rho,double eta, double beta)
+{
+    double c=std::tanh(beta*(rho-eta));
+    double a=std::tanh(eta*beta);
+    double b=std::tanh(beta*(1.0-eta));
+    double rez=-2.0*beta*beta*c*(1.0-c*c)/(a+b);
+    return rez;
+}
+
+
+inline
+double FluidInterpolation(double rho,double q)
+{
+    return q*(1.0-rho)/(q+rho);
+}
+
+inline
+double GradFluidInterpolation(double rho, double q)
+{
+    double tt=q+rho;
+    return -q/tt-q*(1.0-rho)/(tt*tt);
+}
+
+
+}
+
+
+class FilterSolver
+{
+public:
+    FilterSolver(double r_, mfem::ParMesh* pmesh_, int order_=2)
+    {
+        r=r_;
+        order=order_;
+        pmesh=pmesh_;
+        int dim=pmesh->Dimension();
+        sfec=new mfem::H1_FECollection(order, dim);
+        sfes=new mfem::ParFiniteElementSpace(pmesh,sfec,1);
+
+        ifec=new mfem::H1Pos_FECollection(order-1,dim);
+        ifes=new mfem::ParFiniteElementSpace(pmesh,ifec,1);
+
+        dfes=ifes;
+        SetSolver();
+
+
+        K=nullptr;
+        S=nullptr;
+        A=nullptr;
+        pcg=nullptr;
+        prec=nullptr;
+    }
+
+    FilterSolver(double r_, mfem::ParMesh* pmesh_, mfem::ParFiniteElementSpace* dfes_,int order_=2):dfes(dfes_)
+    {
+        r=r_;
+        order=order_;
+        pmesh=pmesh_;
+        int dim=pmesh->Dimension();
+        sfec=new mfem::H1_FECollection(order, dim);
+        sfes=new mfem::ParFiniteElementSpace(pmesh,sfec,1);
+
+        ifec=nullptr;
+        ifes=nullptr;
+
+        SetSolver();
+
+
+        K=nullptr;
+        S=nullptr;
+        A=nullptr;
+        pcg=nullptr;
+        prec=nullptr;
+
+    }
+
+    mfem::ParFiniteElementSpace* GetFilterFES(){return sfes;}
+    mfem::ParFiniteElementSpace* GetDesignFES(){return dfes;}
+
+
+    virtual
+    ~FilterSolver()
+    {
+        delete pcg;
+        delete prec;
+        delete K;
+        delete S;
+        delete A;
+        delete sfes;
+        delete sfec;
+        delete ifes;
+        delete ifec;
+    }
+
+    void Update()
+    {
+        sfes->Update();
+        dfes->Update();
+        AllocSolvers();
+    }
+
+    virtual
+    void Mult(const Vector &x, Vector &y)
+    {
+
+        if(pcg==nullptr){
+            AllocSolvers();
+        }
+
+        //y=bdrc;
+        tmpv.SetSize(y.Size());
+
+        pcg->SetAbsTol(atol);
+        pcg->SetRelTol(rtol);
+        pcg->SetMaxIter(max_iter);
+        pcg->SetPrintLevel(prt_level);
+        S->Mult(x,tmpv);
+        K->EliminateBC(*A,ess_tdofv,bdrc,tmpv);
+        pcg->Mult(tmpv,y);
+
+    }
+
+    virtual
+    void MultTranspose(const Vector &x, Vector &y)
+    {
+        y=0.0;
+        rhsv.SetSize(x.Size()); rhsv=x;
+        tmpv.SetSize(x.Size()); tmpv=0.0;
+        pcg->SetAbsTol(atol);
+        pcg->SetRelTol(rtol);
+        pcg->SetMaxIter(max_iter);
+        pcg->SetPrintLevel(prt_level);
+        K->EliminateBC(*A,ess_tdofv,tmpv,rhsv);
+        pcg->Mult(rhsv,tmpv);
+        S->MultTranspose(tmpv,y);
+    }
+
+    void SetSolver(double rtol_=1e-8, double atol_=1e-12,int miter_=1000, int prt_level_=1)
+    {
+        rtol=rtol_;
+        atol=atol_;
+        max_iter=miter_;
+        prt_level=prt_level_;
+    }
+
+    void AddBC(int id, double val)
+    {
+        bcr[id]=mfem::ConstantCoefficient(val);
+
+        delete pcg;
+        delete prec;
+        delete K;
+        delete S;
+        delete A;
+
+        pcg=nullptr;
+        prec=nullptr;
+        K=nullptr;
+        S=nullptr;
+        A=nullptr;
+
+    }
+
+    void AllocSolvers()
+    {
+        delete pcg;
+        delete prec;
+        delete K;
+        delete S;
+        delete A;
+
+        ess_tdofv.DeleteAll();
+        bdrc.SetSize(sfes->GetTrueVSize()); bdrc=0.0;
+        //set boundary conditions
+        if(bcr.size()!=0)
+        {
+            mfem::ParGridFunction tmpgf(sfes); tmpgf=0.0;
+            for(auto it=bcr.begin();it!=bcr.end();it++)
+            {
+                mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+                ess_bdr=0;
+                ess_bdr[it->first -1]=1;
+                mfem::Array<int> ess_tdof_list;
+                sfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+                ess_tdofv.Append(ess_tdof_list);
+                tmpgf.ProjectBdrCoefficient(it->second,ess_bdr);
+            }
+            tmpgf.GetTrueDofs(bdrc);
+        }
+
+
+        double dr=r/(2.0*sqrt(3.0));
+        mfem::ConstantCoefficient dc(dr*dr);
+
+        mfem::ParBilinearForm* bf=new mfem::ParBilinearForm(sfes);
+        bf->AddDomainIntegrator(new mfem::MassIntegrator());
+        bf->AddDomainIntegrator(new DiffusionIntegrator(dc));
+        bf->Assemble();
+        bf->Finalize();
+        K=bf->ParallelAssemble();
+        delete bf;
+
+        A=K->EliminateRowsCols(ess_tdofv);
+        K->EliminateZeroRows();
+
+        //allocate the CG solver and the preconditioner
+        prec=new mfem::HypreBoomerAMG(*K);
+        pcg=new mfem::CGSolver(pmesh->GetComm());
+        pcg->SetOperator(*K);
+        pcg->SetPreconditioner(*prec);
+
+        mfem::ParMixedBilinearForm* mf=new mfem::ParMixedBilinearForm(dfes,sfes);
+        mf->AddDomainIntegrator(new mfem::MassIntegrator());
+        mf->Assemble();
+        mf->Finalize();
+        S=mf->ParallelAssemble();
+        delete mf;
+
+    }
+
+
+private:
+    mfem::HypreParMatrix* S;
+    mfem::HypreParMatrix* K;
+    mfem::HypreParMatrix* A;
+    mfem::Solver* prec;
+    mfem::CGSolver* pcg;
+
+    mfem::FiniteElementCollection* sfec;
+    mfem::ParFiniteElementSpace* sfes;
+    mfem::FiniteElementCollection* ifec;
+    mfem::ParFiniteElementSpace* ifes;
+
+    mfem::Vector tmpv;
+    mfem::Vector bdrc;//boundary conditions
+    mfem::Vector rhsv;//RHS for the adjoint
+    mfem::Array<int> ess_tdofv; //boundary dofs
+
+    double r;
+    int order;
+
+    mfem::ParMesh* pmesh;
+
+    mfem::ParFiniteElementSpace* dfes;
+
+    std::map<int, mfem::ConstantCoefficient> bcr;
+
+    double atol;
+    double rtol;
+    int max_iter;
+    int prt_level;
+
+};
+
+
+class PVolumeQoIIntegrator:public NonlinearFormIntegrator
+{
+public:
+    PVolumeQoIIntegrator(double eta_, double beta_, int iorder_=4)
+    {
+        eta=eta_;
+        beta=beta_;
+        iorder=iorder_;
+    }
+
+    virtual
+    ~PVolumeQoIIntegrator()
+    {
+
+    }
+
+    virtual
+    double GetElementEnergy(const mfem::FiniteElement &el,
+                            mfem::ElementTransformation &trans,
+                            const mfem::Vector &elfun) override
+    {
+        double energy=0.0;
+        const int ndof = el.GetDof();
+        const int ndim = el.GetDim();
+
+        const mfem::IntegrationRule *ir = NULL;
+        int order = iorder * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::Vector shapef(ndof);
+        double w;
+        double tval;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            tval=shapef*elfun;
+            //trim the density for high-order fields
+            if(tval>1.0){tval=1.0;}
+            else if(tval<0.0){tval=0.0;}
+            w= mfem::PointwiseTrans::HProject(tval,eta,beta);
+            w= ip.weight * trans.Weight() * w;
+            energy = energy + w;
+        }
+        return energy;
+    }
+
+    virtual
+    void AssembleElementVector(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const mfem::Vector & elfun,
+                                       mfem::Vector & elvect) override
+    {
+
+        const int ndof = el.GetDof();
+
+        const mfem::IntegrationRule *ir = NULL;
+        int order = iorder * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        elvect.SetSize(ndof);
+        elvect=0.0;
+
+        mfem::Vector shapef(ndof);
+        double w;
+        double tval;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            tval=shapef*elfun;
+            //trim the density for high-order fields
+            if(tval>1.0){tval=1.0;}
+            else if(tval<0.0){tval=0.0;}
+
+            w= mfem::PointwiseTrans::HGrad(tval,eta,beta);
+            w= ip.weight * trans.Weight() * w;
+            elvect.Add(w,shapef);
+        }
+
+    }
+
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                      mfem::ElementTransformation & trans,
+                                      const mfem::Vector & elfun,
+                                      mfem::DenseMatrix & elmat) override
+     {
+         const int ndof = el.GetDof();
+
+         const mfem::IntegrationRule *ir = NULL;
+         int order = iorder * trans.OrderGrad(&el) - 1; // correct order?
+         ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+         elmat.SetSize(ndof);
+         elmat=0.0;
+
+         mfem::Vector shapef(ndof);
+         double w;
+         double tval;
+         for (int i = 0; i < ir -> GetNPoints(); i++)
+         {
+             const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+             trans.SetIntPoint(&ip);
+             el.CalcShape(ip,shapef);
+             tval=shapef*elfun;
+             //trim the density for high-order fields
+             if(tval>1.0){tval=1.0;}
+             else if(tval<0.0){tval=0.0;}
+             w=mfem::PointwiseTrans::HHess(tval,eta,beta);
+             w= ip.weight * trans.Weight() * w;
+             AddMult_a_VVt(w, shapef, elmat);
+         }
+
+    }
+
+
+private:
+    double eta;
+    double beta;
+    int iorder;
+
+};
+
+
+class VolumeQoIIntegrator:public NonlinearFormIntegrator
+{
+public:
+    VolumeQoIIntegrator(int iorder_=4)
+    {
+        iorder=iorder_;
+    }
+
+    virtual
+    ~VolumeQoIIntegrator()
+    {
+
+    }
+
+    virtual
+    double GetElementEnergy(const mfem::FiniteElement &el,
+                            mfem::ElementTransformation &trans,
+                            const mfem::Vector &elfun) override
+    {
+        double energy=0.0;
+        const int ndof = el.GetDof();
+        const int ndim = el.GetDim();
+
+        const mfem::IntegrationRule *ir = NULL;
+        int order = iorder * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::Vector shapef(ndof);
+        double w;
+        double tval;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            tval=shapef*elfun;
+            //trim the density for high-order fields
+            if(tval>1.0){tval=1.0;}
+            else if(tval<0.0){tval=0.0;}
+            w= tval;
+            w= ip.weight * trans.Weight() * w;
+            energy = energy + w;
+        }
+        return energy;
+    }
+
+    virtual
+    void AssembleElementVector(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const mfem::Vector & elfun,
+                                       mfem::Vector & elvect) override
+    {
+
+        const int ndof = el.GetDof();
+
+        const mfem::IntegrationRule *ir = NULL;
+        int order = iorder * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        elvect.SetSize(ndof);
+        elvect=0.0;
+
+        mfem::Vector shapef(ndof);
+        double w;
+        double tval;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            tval=shapef*elfun;
+            //trim the density for high-order fields
+            if(tval>1.0){tval=1.0;}
+            else if(tval<0.0){tval=0.0;}
+
+            w= 1.0;
+            w= ip.weight * trans.Weight() * w;
+            elvect.Add(w,shapef);
+        }
+
+    }
+
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                      mfem::ElementTransformation & trans,
+                                      const mfem::Vector & elfun,
+                                      mfem::DenseMatrix & elmat) override
+     {
+         const int ndof = el.GetDof();
+
+         const mfem::IntegrationRule *ir = NULL;
+         int order = iorder * trans.OrderGrad(&el) - 1; // correct order?
+         ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+         elmat.SetSize(ndof);
+         elmat=0.0;
+
+         mfem::Vector shapef(ndof);
+         double w;
+         double tval;
+         for (int i = 0; i < ir -> GetNPoints(); i++)
+         {
+             const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+             trans.SetIntPoint(&ip);
+             el.CalcShape(ip,shapef);
+             tval=shapef*elfun;
+             //trim the density for high-order fields
+             if(tval>1.0){tval=1.0;}
+             else if(tval<0.0){tval=0.0;}
+             w=0.0;
+             w= ip.weight * trans.Weight() * w;
+             AddMult_a_VVt(w, shapef, elmat);
+         }
+
+    }
+
+
+private:
+    double eta;
+    double beta;
+    int iorder;
+
+};
+
+
+class VolumeQoI{
+public:
+    VolumeQoI(mfem::ParFiniteElementSpace* fes_)
+    {
+        fes=fes_;
+        nf=nullptr;
+    }
+
+    ~VolumeQoI()
+    {
+        delete nf;
+    }
+
+    void SetProjection(double eta_,double beta_)
+    {
+    }
+
+    void Update()
+    {
+        delete nf;
+        nf=nullptr;
+        fes->Update();
+    }
+
+    /// Input: true design vector
+    double Eval(mfem::Vector& design_)
+    {
+        Alloc();
+        return nf->GetEnergy(design_);
+    }
+
+    void Grad(mfem::Vector& design_,mfem::Vector& grad)
+    {
+        Alloc();
+        nf->Mult(design_,grad);
+    }
+
+private:
+    mfem::ParNonlinearForm* nf;
+    mfem::ParFiniteElementSpace* fes;
+
+    void Alloc()
+    {
+        if(nf==nullptr)
+        {
+            nf=new mfem::ParNonlinearForm(fes);
+            nf->AddDomainIntegrator(new mfem::VolumeQoIIntegrator());
+        }
+    }
+};
+
+//computes the volume of a projected field
+class PVolumeQoI
+{
+public:
+    PVolumeQoI(mfem::ParFiniteElementSpace* fes_)
+    {
+        fes=fes_;
+        nf=nullptr;
+        SetProjection(0.5,8.0);
+    }
+
+    ~PVolumeQoI()
+    {
+        delete nf;
+    }
+
+    void SetProjection(double eta_,double beta_)
+    {
+        eta=eta_;
+        beta=beta_;
+        delete nf;
+        nf=nullptr;
+    }
+
+    void Update()
+    {
+        delete nf;
+        nf=nullptr;
+        fes->Update();
+    }
+
+    /// Input: true design vector
+    double Eval(mfem::Vector& design_)
+    {
+        Alloc();
+        return nf->GetEnergy(design_);
+    }
+
+    void Grad(mfem::Vector& design_,mfem::Vector& grad)
+    {
+        Alloc();
+        nf->Mult(design_,grad);
+    }
+
+private:
+    double eta;
+    double beta;
+    mfem::ParNonlinearForm* nf;
+    mfem::ParFiniteElementSpace* fes;
+
+    void Alloc()
+    {
+        if(nf==nullptr)
+        {
+            nf=new mfem::ParNonlinearForm(fes);
+            nf->AddDomainIntegrator(new mfem::PVolumeQoIIntegrator(eta,beta));
+        }
+    }
+};
+
+}
+
+#endif

miniapps/shifted/nodal_grad_test.cpp

@@ -0,0 +1,156 @@
+#include "mfem.hpp"
+#include "integ_algoim.hpp"
+#include "shape_grad.hpp"
+#include "marking.hpp"
+
+using namespace mfem;
+using namespace std;
+
+//Level set function for sphere in 3D and circle in 2D
+double sphere_ls(const Vector &x)
+{
+   double r2= x*x;
+   return -sqrt(r2)+1.0;//the radius is 1.0
+}
+
+//Level set function for a sinusoidal wave.
+//Resulting zero isocontour is at y=0.5-(0.1*sin(3*pi*x+pi/2))
+double sinusoidal_ls(const Vector &x)
+{
+   double a1 = 20., a2 = 2., a3 = 3.;
+   return tanh(a1*(x(1)-0.5) + a2*sin(a3*(x(0)-0.5)*M_PI));
+}
+
+
+int main(int argc, char *argv[])
+{
+   //   Parse command-line options
+   const char *mesh_file = "../../data/star-q3.mesh";
+   int ser_ref_levels = 1;
+   int order = 1;
+   bool visualization = true;
+   int print_level = 0;
+   int ls_type = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ls_type,
+                  "-ls",
+                  "--ls-type",
+                  "Level set type: 1: circle, 2 sinusoidal wave");
+   args.AddOption((&print_level), "-prt", "--print-level", "Print level.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(std::cout);
+      return 1;
+   }
+   args.PrintOptions(std::cout);
+
+   //   Read the (serial) mesh from the given mesh file.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+   //    Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   mesh->EnsureNodes();
+   GridFunction* nodes=mesh->GetNodes();
+
+   // Define the finite element space for the level-set function.
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(mesh, &fec, 1, Ordering::byVDIM);
+   int glob_size = fespace.GetTrueVSize();
+   std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+
+   // Define the level set grid function
+   GridFunction x(&fespace);
+
+   // Define the level set coefficient
+   Coefficient *ls_coeff = nullptr;
+   if (ls_type == 1)
+   {
+      ls_coeff=new FunctionCoefficient(sphere_ls);
+   }
+   else if (ls_type == 2)
+   {
+      ls_coeff=new FunctionCoefficient(sinusoidal_ls);
+   }
+   else
+   {
+      MFEM_ABORT("Level set coefficient not defined");
+   }
+
+   // Project the coefficient onto the LS grid function
+   x.ProjectCoefficient(*ls_coeff);
+
+   ConstantCoefficient one(1.0);
+   Array<int> el_marks; el_marks.SetSize(mesh->GetNE());
+   for(int i=0;i<mesh->GetNE();i++){ el_marks[i]= ElementMarker::SBElementType::CUT;}
+
+
+   NonlinearForm nf(&fespace);
+   nf.AddDomainIntegrator(new DVolShapeIntegrator(one,el_marks));
+
+   double vol=nf.GetEnergy(x);
+   Vector grad(x.Size());
+   nf.Mult(x,grad);
+   std::cout<<setprecision(8);
+   grad.Print(std::cout,x.Size());
+
+   std::cout<<"FD+"<<std::endl;
+   double dd=0.00000001;
+   for(int i=0;i<x.Size();i++){
+       x[i]=x[i]-dd;
+       double vol1=nf.GetEnergy(x);
+       x[i]=x[i]+dd;
+       std::cout<<(vol1-vol)/dd<<" ";
+   }
+   std::cout<<std::endl;
+
+   std::cout<<"vol="<<vol<<std::endl;
+
+   Vector di(nodes->Size()); di=0.0;
+   Vector si(nodes->Size()); si=0.0;
+
+   for(int i=0;i<nodes->Size();i++)
+   {
+       (*nodes)[i]=(*nodes)[i]+dd;
+       mesh->DeleteGeometricFactors();
+       double vol1=nf.GetEnergy(x);
+       si[i]=(vol1-vol)/dd;
+       (*nodes)[i]=(*nodes)[i]-dd;
+   }
+
+   std::cout<<"si"<<std::endl;
+   si.Print(std::cout,2);
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaViewDistance", mesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.RegisterField("x", &x);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+
+
+
+   delete ls_coeff;
+   delete mesh;
+
+   return 0;
+}

miniapps/shifted/shape_grad.cpp

@@ -0,0 +1,627 @@
+#include "shape_grad.hpp"
+#include "marking.hpp"
+
+
+#ifdef MFEM_USE_ALGOIM
+#include "integ_algoim.hpp"
+#endif
+
+namespace mfem{
+
+#ifdef MFEM_USE_ALGOIM
+
+void DVolShapeIntegrator::AssembleElementVector(const FiniteElement &el,
+                                          ElementTransformation &Tr,
+                                          const Vector &elfun, Vector &elvect)
+{
+
+    if((*elem_markers)[Tr.ElementNo]==ElementMarker::SBElementType::CUT)
+    {
+        int ndof = el.GetDof();
+        int ndim = Tr.GetSpaceDim();
+
+        elvect.SetSize(ndof); elvect=0.0;
+
+        DenseMatrix bmat(ndof,ndim); //gradients of the shape functions in isoparametric space
+        Vector tnormal(ndim); //normal to the level set in true space
+        Vector shf(ndof);
+        DenseMatrix gradn(ndof,ndim); //nodal gradient values
+
+        DenseMatrix proj(ndim*ndof,ndof);
+        el.ProjectGrad(el,Tr,proj);
+        DenseMatrix ngrad[ndim];
+        for(int i=0;i<ndim;i++){
+            ngrad[i].SetSize(ndof,ndof);
+            proj.GetSubMatrix(i*ndof,(i+1)*ndof,0,ndof,ngrad[i]);
+        }
+
+        {
+            Vector gradv; gradv.SetDataAndSize(gradn.GetData(),ndim*ndof);
+            proj.Mult(elfun,gradv);
+        }
+
+        int order;
+        if(lorder>0){
+            order=lorder;
+        }else{
+            order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+        }
+
+        order=10;
+
+        AlgoimIntegrationRule air(order, el, Tr, elfun);
+        const IntegrationRule *ir = air.GetVolumeIntegrationRule();
+
+        double w;
+        double f;
+        DenseMatrix pmat(ndof,ndim);
+
+        //evaluate the gradients with respect to nodal displacements
+        /*
+        DenseMatrix gp(ndof,ndim); gp=0.0;
+        {
+            for (int j = 0; j < ir->GetNPoints(); j++)
+            {
+                const IntegrationPoint &ip = ir->IntPoint(j);
+                Tr.SetIntPoint(&ip);
+                el.CalcDShape(ip,bmat);
+                Mult(bmat,Tr.AdjugateJacobian(),pmat);
+                w = ip.weight ;
+                f = coeff->Eval(Tr,ip);
+                gp.Add(w*f,pmat);
+            }
+            std::cout<<"GP="<<std::endl;
+            gp.PrintMatlab(std::cout);
+        }*/
+
+        Vector gradf(ndim); gradf=0.0;
+        Vector lap(ndof);
+        Vector tv(ndof);
+        Vector gv;
+        for (int j = 0; j < ir->GetNPoints(); j++)
+        {
+           const IntegrationPoint &ip = ir->IntPoint(j);
+           Tr.SetIntPoint(&ip);
+           el.CalcPhysDShape(Tr,bmat);
+
+           f = coeff->Eval(Tr,ip);
+           w = ip.weight*Tr.Weight();
+
+
+           for(int d=0;d<ndim;d++){
+               gv.SetDataAndSize(bmat.GetData()+d*ndof,ndof);
+               ngrad[d].MultTranspose(gv,tv);
+               elvect.Add(w*f,tv);
+           }
+
+           bmat.Mult(gradf,tv);
+           elvect.Add(w,tv);
+        }
+
+
+        elvect.Print(std::cout);
+
+    }
+    else
+    {
+        elvect.SetSize(elfun.Size());
+        elvect=0.0;
+    }
+
+
+}
+
+double DVolShapeIntegrator::GetElementEnergy(const FiniteElement &el,
+                                            ElementTransformation &Tr,
+                                            const Vector &elfun)
+{
+    if((*elem_markers)[Tr.ElementNo]==ElementMarker::SBElementType::INSIDE){
+        double val=0.0;
+        const IntegrationRule * ir = nullptr;
+        int order;
+        if(lorder>0){
+            order=lorder;
+        }else{
+            order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+        }
+        ir = &IntRules.Get(el.GetGeomType(), order);
+        double w;
+        double f;
+        for(int i=0; i < ir->GetNPoints(); i++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w = Tr.Weight();
+            w = ip.weight * w;
+            f = coeff->Eval(Tr,ip);
+            val = val + f*w;
+
+        }
+
+        return val;
+    }else
+    if((*elem_markers)[Tr.ElementNo]==ElementMarker::SBElementType::OUTSIDE)
+    {
+        return 0.0;
+    }
+
+    // integrate cut element
+    double val = 0.0;
+    int order;
+    if(lorder>0){
+        order=lorder;
+    }else{
+        order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+    }
+    AlgoimIntegrationRule air(order, el, Tr, elfun);
+    const IntegrationRule * ir = air.GetVolumeIntegrationRule();
+    for (int j = 0; j < ir->GetNPoints(); j++)
+    {
+       const IntegrationPoint &ip = ir->IntPoint(j);
+       Tr.SetIntPoint(&ip);
+       val += ip.weight * Tr.Weight();
+    }
+
+    return val;
+}
+
+void VolShapeIntegrator::AssembleElementVector(const FiniteElement &el,
+                                          ElementTransformation &Tr,
+                                          const Vector &elfun, Vector &elvect)
+{
+
+    if((*elem_markers)[Tr.ElementNo]==ElementMarker::SBElementType::CUT)
+    {
+        int ndof = el.GetDof();
+        int ndim = Tr.GetSpaceDim();
+
+        elvect.SetSize(ndof); elvect=0.0;
+
+        DenseMatrix bmat(ndof,ndim); //gradients of the shape functions in isoparametric space
+        DenseMatrix pmat(ndof,ndim);
+        Vector inormal(ndim); //normal to the level set in isoparametric space
+        Vector tnormal(ndim);
+        Vector shf(ndof);
+
+        int order;
+        if(lorder>0){
+            order=lorder;
+        }else{
+            order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+        }
+
+        AlgoimIntegrationRule air(order, el, Tr, elfun);
+        const IntegrationRule * ir = air.GetSurfaceIntegrationRule();
+
+        double w;
+        double f;
+
+        for (int j = 0; j < ir->GetNPoints(); j++)
+        {
+           const IntegrationPoint &ip = ir->IntPoint(j);
+           Tr.SetIntPoint(&ip);
+           el.CalcDShape(ip,bmat);
+           Mult(bmat, Tr.AdjugateJacobian(), pmat);
+
+
+           //compute the normal to the LS in isoparametric space
+           bmat.MultTranspose(elfun,inormal);
+           pmat.MultTranspose(elfun,tnormal);
+           std::cout<<"tn="<<tnormal.Norml2()<<" in="<<inormal.Norml2()<<" ww"<< tnormal.Norml2()/inormal.Norml2()<<std::endl;
+           std::cout<<"sqrt(w)="<<sqrt(Tr.Weight())<<" w^2/3="<<pow(Tr.Weight(), 2.0/3.0)<<std::endl;
+
+           el.CalcPhysDShape(Tr,pmat);
+           pmat.MultTranspose(elfun,tnormal);
+           std::cout<<"tn="<<tnormal.Norml2()<<std::endl;
+
+           //w=ip.weight*(tnormal.Norml2())/inormal.Norml2();
+           //w=ip.weight/tnormal.Norml2();
+           if(ndim==2) { w = ip.weight * sqrt(Tr.Weight())/tnormal.Norml2();} /// inormal.Norml2();
+           else { w= ip.weight * pow(Tr.Weight(), 2.0/3.0)/tnormal.Norml2();}
+           el.CalcPhysShape(Tr,shf);
+           f = coeff->Eval(Tr,ip);
+           elvect.Add(w * f , shf);
+        }
+    }
+    else
+    {
+        elvect.SetSize(elfun.Size());
+        elvect=0.0;
+    }
+
+
+}
+
+
+double VolShapeIntegrator::GetElementEnergy(const FiniteElement &el,
+                                            ElementTransformation &Tr,
+                                            const Vector &elfun)
+{
+    if((*elem_markers)[Tr.ElementNo]==ElementMarker::SBElementType::INSIDE){
+        double val=0.0;
+        const IntegrationRule * ir = nullptr;
+        int order;
+        if(lorder>0){
+            order=lorder;
+        }else{
+            order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+        }
+        ir = &IntRules.Get(el.GetGeomType(), order);
+        double w;
+        double f;
+        for(int i=0; i < ir->GetNPoints(); i++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w = Tr.Weight();
+            w = ip.weight * w;
+            f = coeff->Eval(Tr,ip);
+            val = val + f*w;
+
+        }
+
+        return val;
+    }else
+    if((*elem_markers)[Tr.ElementNo]==ElementMarker::SBElementType::OUTSIDE)
+    {
+        return 0.0;
+    }
+
+    // integrate cut element
+    double val = 0.0;
+    int order;
+    if(lorder>0){
+        order=lorder;
+    }else{
+        order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+    }
+    AlgoimIntegrationRule air(order, el, Tr, elfun);
+    const IntegrationRule * ir = air.GetVolumeIntegrationRule();
+    for (int j = 0; j < ir->GetNPoints(); j++)
+    {
+       const IntegrationPoint &ip = ir->IntPoint(j);
+       Tr.SetIntPoint(&ip);
+       val += ip.weight * Tr.Weight();
+    }
+
+    return val;
+}
+
+double SurfShapeIntegrator::GetElementEnergy(const FiniteElement &el,
+                                            ElementTransformation &Tr,
+                                            const Vector &elfun)
+{
+    if((*elem_markers)[Tr.ElementNo]!=ElementMarker::SBElementType::CUT)
+    {
+        return 0.0;
+    }
+
+    // integrate cut element
+    double val = 0.0;
+    int order;
+    if(lorder>0){
+        order=lorder;
+    }else{
+        order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+    }
+
+    AlgoimIntegrationRule air(order, el, Tr, elfun);
+    const IntegrationRule * ir = air.GetSurfaceIntegrationRule();
+
+    int ndof = el.GetDof();
+    int ndim = Tr.GetSpaceDim();
+
+    DenseMatrix bmat(ndof,ndim); //gradients of the shape functions in isoparametric space
+    DenseMatrix pmat(ndof,ndim); //gradients of the shape functions in physical space
+    Vector inormal(ndim); //normal to the level set in isoparametric space
+    Vector tnormal(ndim);
+
+    double w;
+    double f;
+
+    for (int j = 0; j < ir->GetNPoints(); j++)
+    {
+       const IntegrationPoint &ip = ir->IntPoint(j);
+       Tr.SetIntPoint(&ip);
+       el.CalcDShape(ip,bmat);
+       Mult(bmat, Tr.AdjugateJacobian(), pmat);
+       //compute the normal to the LS in isoparametric space
+       bmat.MultTranspose(elfun,inormal);
+       //compute the normal to the LS in physical space
+       pmat.MultTranspose(elfun,tnormal);
+
+       //std::cout<<" sca1="<< tnormal.Norml2()/ inormal.Norml2()<< " (adjJ)="<<Tr.AdjugateJacobian().Det() <<" |J|=" <<Tr.Weight()<<" det(J)="<< Tr.Jacobian().Det()<<std::endl;
+
+       w = ip.weight * tnormal.Norml2()/ inormal.Norml2();
+       f = coeff->Eval(Tr,ip);
+
+       val=val + w*f;
+    }
+
+    return val;
+}
+
+double SurfMeanCurv3D(Vector& sh, Vector& gradv, Vector& dgradx, Vector& dgrady, Vector& dgradz)
+{
+    int ndof=sh.Size();
+    Vector tv;
+
+    tv.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+    double tx=sh*tv;
+    tv.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+    double ty=sh*tv;
+    tv.SetDataAndSize(gradv.GetData()+2*ndof,ndof);
+    double tz=sh*tv;
+
+    tv.SetDataAndSize(dgradx.GetData()+0*ndof,ndof);
+    double txx=sh*tv;
+    tv.SetDataAndSize(dgradx.GetData()+1*ndof,ndof);
+    double txy=sh*tv;
+    tv.SetDataAndSize(dgradx.GetData()+2*ndof,ndof);
+    double txz=sh*tv;
+
+    tv.SetDataAndSize(dgrady.GetData()+1*ndof,ndof);
+    double tyy=sh*tv;
+    tv.SetDataAndSize(dgrady.GetData()+2*ndof,ndof);
+    double tyz=sh*tv;
+
+    tv.SetDataAndSize(dgradz.GetData()+2*ndof,ndof);
+    double tzz=sh*tv;
+
+    double nr=sqrt(tx*tx+ty*ty+tz*tz);
+    double rez=tx*tx*(tyy+tzz)+ty*ty*(txx+tzz)+tz*tz*(txx+tyy);
+    rez=rez-2.0*tx*ty*txy-2.0*tx*tz*txz-2.0*ty*tz*tyz;
+    rez=rez/(nr*nr*nr);
+    return rez;
+}
+
+double SurfMeanCurv2D(Vector& sh, Vector& gradv, Vector& dgradx, Vector& dgrady, Vector& dgradz)
+{
+    int ndof=sh.Size();
+    Vector tv;
+
+    tv.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+    double tx=sh*tv;
+    tv.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+    double ty=sh*tv;
+
+    tv.SetDataAndSize(dgradx.GetData()+0*ndof,ndof);
+    double txx=sh*tv;
+    tv.SetDataAndSize(dgradx.GetData()+1*ndof,ndof);
+    double txy=sh*tv;
+
+    tv.SetDataAndSize(dgrady.GetData()+1*ndof,ndof);
+    double tyy=sh*tv;
+
+
+    double nr=sqrt(tx*tx+ty*ty);
+    double rez=tx*tx*(tyy)+ty*ty*(txx);
+    rez=rez-2.0*tx*ty*txy;
+    rez=rez/(nr*nr*nr);
+    return rez;
+}
+
+void SurfShapeIntegrator::AssembleElementVector(const FiniteElement &el,
+                                          ElementTransformation &Tr,
+                                          const Vector &elfun, Vector &elvect)
+{
+    elvect.SetSize(elfun.Size()); elvect=0.0;
+    if((*elem_markers)[Tr.ElementNo]!=ElementMarker::SBElementType::CUT)
+    {
+        return;
+    }
+
+    int ndof = el.GetDof();
+    int ndim = Tr.GetSpaceDim();
+
+    DenseMatrix proj(ndim*ndof,ndof);
+    el.ProjectGrad(el,Tr,proj);
+
+    Vector gradv(ndim*ndof);
+    proj.Mult(elfun,gradv);
+
+    Vector gradx; gradx.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+    Vector grady; grady.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+    Vector gradz;
+    if(ndim==3){
+        gradz.SetDataAndSize(gradv.GetData()+2*ndof,ndof);
+    }else{
+        gradz.SetSize(ndof); gradz=0.0;
+    }
+
+    Vector dgradx(ndim*ndof); proj.Mult(gradx,dgradx);
+    Vector dgrady(ndim*ndof); proj.Mult(grady,dgrady);
+    Vector dgradz(ndim*ndof); proj.Mult(gradz,dgradz);
+
+    DenseMatrix bmat(ndof,ndim); //gradients of the shape functions in isoparametric space
+    Vector inormal(ndim); //normal to the level set in isoparametric space
+    Vector shf(ndof);
+
+    int order;
+    if(lorder>0){
+        order=lorder;
+    }else{
+        order = 2 * el.GetOrder() +Tr.OrderGrad(&el);
+    }
+
+    AlgoimIntegrationRule air(order, el, Tr, elfun);
+    const IntegrationRule * ir = air.GetSurfaceIntegrationRule();
+
+    double w;
+    double f;
+    Vector H;
+    Vector nn(ndim);
+    Vector gr(ndim); gr=0.0;
+
+    //evaluate the curvature at the integrations points
+    MeanCurvImplicitFunction(elfun,el,Tr,*ir,H);
+
+    for (int j = 0; j < ir->GetNPoints(); j++)
+    {
+       const IntegrationPoint &ip = ir->IntPoint(j);
+       Tr.SetIntPoint(&ip);
+       el.CalcDShape(ip,bmat);
+       //compute the normal to the LS in isoparametric space
+       bmat.MultTranspose(elfun,inormal);
+       w = ip.weight * Tr.Weight() / inormal.Norml2();
+
+       el.CalcPhysDShape(Tr,bmat);
+       bmat.MultTranspose(elfun,inormal);
+       if(ndim==2) { w = ip.weight * sqrt(Tr.Weight())/inormal.Norml2();} /// inormal.Norml2();
+       else { w= ip.weight * pow(Tr.Weight(), 2.0/3.0)/inormal.Norml2();}
+
+
+
+       el.CalcPhysShape(Tr,shf);
+
+       double H1;
+       //Compute the mean curvature H
+       if(ndim==3){
+           H1=SurfMeanCurv3D(shf,gradv,dgradx,dgrady,dgradz);
+       }else{
+           H1=SurfMeanCurv2D(shf,gradv,dgradx,dgrady,dgradz);
+       }
+
+       nn[0]=shf*gradx;
+       nn[1]=shf*grady;
+       if(ndim==3){ nn[2]=shf*gradz;}
+
+       std::cout<<"nn="<<nn.Norml2()<<" H="<<H[j]<<" H1="<<H1<<std::endl;
+
+       f = coeff->Eval(Tr,ip);
+
+       if(gradco!=nullptr){
+           gradco->Eval(gr,Tr,ip);
+       }
+
+       gr.Print(std::cout); std::cout<<" dot="<<nn*gr/nn.Norml2()<<std::endl;
+
+       //elvect.Add(w*(-H*f),shf);
+       elvect.Add(w*(-H[j]*f-nn*gr/nn.Norml2()),shf);
+
+
+
+    }
+}
+
+
+/// Evaluates the mean curvature div(n/|n|) for all integration points
+/// elfun - nodal values of a level-set function
+void MeanCurvImplicitFunction(const Vector& elfun,
+                    const FiniteElement & el,
+                    ElementTransformation &T,
+                    const IntegrationRule& ir,
+                    Vector& H)
+{
+    H.SetSize(ir.GetNPoints()); H=0.0;
+
+    int ndim = el.GetDim();
+    int ndof = el.GetDof();
+
+    DenseMatrix proj(ndim*ndof,ndof);
+    el.ProjectGrad(el,T,proj);
+
+    Vector gradv(ndim*ndof);
+    proj.Mult(elfun,gradv);
+
+    Vector grad[ndim];
+    for(int i=0;i<ndim;i++){
+        grad[i].SetDataAndSize(gradv.GetData()+i*ndof,ndof);
+    }
+
+    Vector dgrad[ndim];
+    for(int i=0;i<ndim;i++){
+        dgrad[i].SetSize(ndim*ndof);
+        proj.Mult(grad[i],dgrad[i]);
+    }
+
+    Vector sh(ndof);
+    Vector tv;
+
+    if(ndim==2){
+    for (int j = 0; j < ir.GetNPoints(); j++)
+    {
+        const IntegrationPoint &ip = ir.IntPoint(j);
+        T.SetIntPoint(&ip);
+        el.CalcPhysShape(T,sh);
+
+        tv.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+        double tx=sh*tv;
+        tv.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+        double ty=sh*tv;
+
+        tv.SetDataAndSize(dgrad[0].GetData()+0*ndof,ndof);
+        double txx=sh*tv;
+        tv.SetDataAndSize(dgrad[0].GetData()+1*ndof,ndof);
+        double txy=sh*tv;
+
+        tv.SetDataAndSize(dgrad[1].GetData()+1*ndof,ndof);
+        double tyy=sh*tv;
+
+        double nr=sqrt(tx*tx+ty*ty);
+        double rez=tx*tx*(tyy)+ty*ty*(txx);
+        rez=rez-2.0*tx*ty*txy;
+        rez=rez/(nr*nr*nr);
+        H[j]=rez;
+
+    }}else{//ndim=3
+    for (int j = 0; j < ir.GetNPoints(); j++){
+        const IntegrationPoint &ip = ir.IntPoint(j);
+        T.SetIntPoint(&ip);
+        el.CalcPhysShape(T,sh);
+
+        tv.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+        double tx=sh*tv;
+        tv.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+        double ty=sh*tv;
+        tv.SetDataAndSize(gradv.GetData()+2*ndof,ndof);
+        double tz=sh*tv;
+
+        tv.SetDataAndSize(dgrad[0].GetData()+0*ndof,ndof);
+        double txx=sh*tv;
+        tv.SetDataAndSize(dgrad[0].GetData()+1*ndof,ndof);
+        double txy=sh*tv;
+        tv.SetDataAndSize(dgrad[0].GetData()+2*ndof,ndof);
+        double txz=sh*tv;
+
+        tv.SetDataAndSize(dgrad[1].GetData()+1*ndof,ndof);
+        double tyy=sh*tv;
+        tv.SetDataAndSize(dgrad[1].GetData()+2*ndof,ndof);
+        double tyz=sh*tv;
+
+        tv.SetDataAndSize(dgrad[2].GetData()+2*ndof,ndof);
+        double tzz=sh*tv;
+
+        double nr=sqrt(tx*tx+ty*ty+tz*tz);
+        double rez=tx*tx*(tyy+tzz)+ty*ty*(txx+tzz)+tz*tz*(txx+tyy);
+        rez=rez-2.0*tx*ty*txy-2.0*tx*tz*txz-2.0*ty*tz*tyz;
+        rez=rez/(nr*nr*nr);
+        H[j]=rez;
+    }}
+}
+
+/// Evaluates the mean curvature div(n/|n|) for all integration points
+/// elfun - nodal values of a level-set function
+void MeanCurvImplicitFunction(int elno,
+                              GridFunction& gf,
+                              const IntegrationRule& ir,
+                              Vector& H)
+{
+    const FiniteElement* el= gf.FESpace()->GetFE(elno);
+    int ndof = el->GetDof();
+
+    Vector elfun(ndof);
+    Array<int> dofs;
+    gf.FESpace()->GetElementDofs(elno, dofs);
+    gf.FESpace()->DofsToVDofs(dofs);
+    gf.GetSubVector(dofs,elfun);
+
+    ElementTransformation* T=gf.FESpace()->GetElementTransformation(elno);
+    MeanCurvImplicitFunction(elfun,*el,*T,ir,H);
+}
+
+
+
+#endif
+
+}

miniapps/shifted/shape_grad.hpp

@@ -0,0 +1,291 @@
+#ifndef SHAPE_GRAD_HPP
+#define SHAPE_GRAD_HPP
+
+
+#include "mfem.hpp"
+
+#ifdef MFEM_USE_ALGOIM
+namespace mfem{
+
+/// Volumetric shape integrator - discrete
+class DVolShapeIntegrator: public NonlinearFormIntegrator
+{
+
+public:
+    DVolShapeIntegrator(Coefficient& coeff_, Array<int> &elem_markers_, int order_=-1)
+    {
+        coeff=&coeff_;
+        elem_markers=&elem_markers_;
+        lorder=order_;
+    }
+
+    virtual
+    ~DVolShapeIntegrator()
+    {
+
+    }
+
+    /// Perform the local action of the NonlinearFormIntegrator
+    virtual void AssembleElementVector(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       const Vector &elfun, Vector &elvect);
+
+    /// Compute the local energy
+    virtual double GetElementEnergy(const FiniteElement &el,
+                                    ElementTransformation &Tr,
+                                    const Vector &elfun);
+
+private:
+    Coefficient* coeff;
+    Array<int>* elem_markers;
+    int lorder;
+};
+
+
+/// Volumetric shape integrator
+class VolShapeIntegrator: public NonlinearFormIntegrator
+{
+
+public:
+    VolShapeIntegrator(Coefficient& coeff_, Array<int> &elem_markers_, int order_=-1)
+    {
+        coeff=&coeff_;
+        elem_markers=&elem_markers_;
+        lorder=order_;
+    }
+
+    virtual
+    ~VolShapeIntegrator()
+    {
+
+    }
+
+    /// Perform the local action of the NonlinearFormIntegrator
+    virtual void AssembleElementVector(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       const Vector &elfun, Vector &elvect);
+
+    /// Compute the local energy
+    virtual double GetElementEnergy(const FiniteElement &el,
+                                    ElementTransformation &Tr,
+                                    const Vector &elfun);
+
+private:
+    Coefficient* coeff;
+    Array<int>* elem_markers;
+    int lorder;
+};
+
+class VolObjectiveCutA
+{
+public:
+    VolObjectiveCutA():volw(1.0)
+    {
+        volc=&volw;
+    }
+
+    void SetWeight(Coefficient* coeff){
+        volc=coeff;
+    }
+
+    void SetCutIntegrationRules(Array<int>& el_markers_)
+    {
+        marks=&el_markers_;
+    }
+
+
+    double Eval(ParGridFunction& lsf){
+        VolShapeIntegrator* itg=new VolShapeIntegrator(*volc,*marks);
+        ParNonlinearForm* nf=new ParNonlinearForm(lsf.ParFESpace());
+        nf->AddDomainIntegrator(itg);
+        double vol=nf->GetEnergy(lsf.GetTrueVector());
+        delete nf;
+        return vol;
+    }
+
+    void Grad(ParGridFunction& lsf, Vector& grad){
+        VolShapeIntegrator* itg=new VolShapeIntegrator(*volc,*marks);
+        ParNonlinearForm* nf=new ParNonlinearForm(lsf.ParFESpace());
+        nf->AddDomainIntegrator(itg);
+        grad.SetSize(lsf.GetTrueVector().Size());
+        nf->Mult(lsf.GetTrueVector(),grad);
+        delete nf;
+    }
+
+
+private:
+    ConstantCoefficient volw; //volume weight
+    Coefficient* volc; //points either to volw or to user supplied coefficient
+
+    Array<int>* marks;
+
+};
+
+
+
+/// Surface shape integrator
+class SurfShapeIntegrator: public NonlinearFormIntegrator
+{
+public:
+    SurfShapeIntegrator(Coefficient& coeff_, Array<int> &elem_markers_, int order_=-1)
+    {
+        coeff=&coeff_;
+        elem_markers=&elem_markers_;
+        lorder=order_;
+        gradco=nullptr;
+    }
+
+    SurfShapeIntegrator(Coefficient& coeff_, VectorCoefficient& gradco_, Array<int> &elem_markers_, int order_=-1)
+    {
+        coeff=&coeff_;
+        elem_markers=&elem_markers_;
+        lorder=order_;
+        gradco=&gradco_;
+    }
+
+    virtual
+    ~SurfShapeIntegrator(){}
+
+    /// Perform the local action of the NonlinearFormIntegrator
+    virtual void AssembleElementVector(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       const Vector &elfun, Vector &elvect);
+
+    /// Compute the local energy
+    virtual double GetElementEnergy(const FiniteElement &el,
+                                    ElementTransformation &Tr,
+                                    const Vector &elfun);
+
+private:
+    Coefficient* coeff;
+    VectorCoefficient* gradco;
+    Array<int>* elem_markers;
+    int lorder;
+};
+
+
+/// Evaluates div(n/|n|) for a level set function
+class MeanCurvImplicitSurf:public Coefficient
+{
+public:
+    MeanCurvImplicitSurf(GridFunction & gf_)
+    {
+        gf=&gf_;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+    {
+        T.SetIntPoint(&ip);
+
+        const FiniteElement* el= gf->FESpace()->GetFE(T.ElementNo);
+        int ndim = el->GetDim();
+        int ndof = el->GetDof();
+
+        Vector elfun(ndof);
+        Array<int> dofs;
+        gf->FESpace()->GetElementDofs(T.ElementNo, dofs);
+        gf->FESpace()->DofsToVDofs(dofs);
+        gf->GetSubVector(dofs,elfun);
+
+        DenseMatrix proj(ndim*ndof,ndof);
+        el->ProjectGrad(*el,T,proj);
+
+
+        Vector gradv(ndim*ndof);
+        proj.Mult(elfun,gradv);
+        Vector gradx; gradx.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+        Vector grady; grady.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+        Vector gradz;
+        if(ndim==3){gradz.SetDataAndSize(gradv.GetData()+2*ndof,ndof);}
+
+
+        Vector dgradx(ndim*ndof); proj.Mult(gradx,dgradx);
+        Vector dgrady(ndim*ndof); proj.Mult(grady,dgrady);
+        Vector dgradz(ndim*ndof);
+        if(ndim==3){ proj.Mult(gradz,dgradz);}
+
+        Vector sh(ndof);
+        el->CalcPhysShape(T,sh);
+
+        if(ndim==2){
+            Vector tv;
+
+            tv.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+            double tx=sh*tv;
+            tv.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+            double ty=sh*tv;
+
+            tv.SetDataAndSize(dgradx.GetData()+0*ndof,ndof);
+            double txx=sh*tv;
+            tv.SetDataAndSize(dgradx.GetData()+1*ndof,ndof);
+            double txy=sh*tv;
+
+            tv.SetDataAndSize(dgrady.GetData()+1*ndof,ndof);
+            double tyy=sh*tv;
+
+
+            double nr=sqrt(tx*tx+ty*ty);
+            double rez=tx*tx*(tyy)+ty*ty*(txx);
+            rez=rez-2.0*tx*ty*txy;
+            rez=rez/(nr*nr*nr);
+            return rez;
+
+        }else{
+            //ndim==3
+            Vector tv;
+
+            tv.SetDataAndSize(gradv.GetData()+0*ndof,ndof);
+            double tx=sh*tv;
+            tv.SetDataAndSize(gradv.GetData()+1*ndof,ndof);
+            double ty=sh*tv;
+            tv.SetDataAndSize(gradv.GetData()+2*ndof,ndof);
+            double tz=sh*tv;
+
+            tv.SetDataAndSize(dgradx.GetData()+0*ndof,ndof);
+            double txx=sh*tv;
+            tv.SetDataAndSize(dgradx.GetData()+1*ndof,ndof);
+            double txy=sh*tv;
+            tv.SetDataAndSize(dgradx.GetData()+2*ndof,ndof);
+            double txz=sh*tv;
+
+            tv.SetDataAndSize(dgrady.GetData()+1*ndof,ndof);
+            double tyy=sh*tv;
+            tv.SetDataAndSize(dgrady.GetData()+2*ndof,ndof);
+            double tyz=sh*tv;
+
+            tv.SetDataAndSize(dgradz.GetData()+2*ndof,ndof);
+            double tzz=sh*tv;
+
+            double nr=sqrt(tx*tx+ty*ty+tz*tz);
+            double rez=tx*tx*(tyy+tzz)+ty*ty*(txx+tzz)+tz*tz*(txx+tyy);
+            rez=rez-2.0*tx*ty*txy-2.0*tx*tz*txz-2.0*ty*tz*tyz;
+            rez=rez/(nr*nr*nr);
+            return rez;
+        }
+    }
+
+private:
+    GridFunction *gf;
+};
+
+/// Evaluates the mean curvature div(n/|n|) for all integration points
+/// elfun - nodal values of a level-set function
+void MeanCurvImplicitFunction(const Vector& elfun,
+                    const FiniteElement & el,
+                    ElementTransformation &T,
+                    const IntegrationRule& ir,
+                    Vector& H);
+
+/// Evaluates the mean curvature div(n/|n|) for all integration points
+/// elfun - nodal values of a level-set function
+void MeanCurvImplicitFunction(int elno,
+                              GridFunction& gf,
+                              const IntegrationRule& ir,
+                              Vector& H);
+
+}
+
+#endif
+
+#endif

miniapps/shifted/shape_grad_test.cpp

@@ -0,0 +1,242 @@
+#include <fstream>
+#include <iostream>
+#include "../common/mfem-common.hpp"
+#include "marking.hpp"
+
+#ifdef MFEM_USE_ALGOIM
+#include "shape_grad.hpp"
+#endif
+
+using namespace mfem;
+using namespace std;
+
+
+const double radius = 0.4;
+
+double sphere_ls(const Vector &x)
+{
+   const int dim = x.Size();
+   const double xc = x(0) - 0.5;
+   const double yc = (dim > 1) ? x(1) - 0.5 : 0.0;
+   const double zc = (dim > 2) ? x(2) - 0.5 : 0.0;
+   const double r = sqrt(xc*xc + yc*yc + zc*zc);
+
+   return (r >= radius) ? -1.0 : 1.0;
+}
+
+double exact_dist_sphere(const Vector &x)
+{
+   const int dim = x.Size();
+   const double xc = x(0) - 0.5;
+   const double yc = (dim > 1) ? x(1) - 0.5 : 0.0;
+   const double zc = (dim > 2) ? x(2) - 0.5 : 0.0;
+   const double r = sqrt(xc*xc + yc*yc + zc*zc);
+
+   //return 3*(radius*radius*radius-r*r*r);
+   return 3*(radius*radius*radius-r*r*r);
+}
+
+double Gyroid(const Vector &xx)
+{
+   const double period = 2.0 * M_PI;
+   double x = xx[0]*period;
+   double y = xx[1]*period;
+   double z = (xx.Size()==3) ? xx[2]*period : 0.0;
+
+   return std::sin(x)*std::cos(y) +
+          std::sin(y)*std::cos(z) +
+          std::sin(z)*std::cos(x) + x*x+y*y;
+}
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int solver_type = 0;
+   int problem = 1;
+   int rs_levels = 2;
+   int order = 2;
+   double t_param = 1.0;
+   const char *device_config = "cpu";
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&solver_type, "-s", "--solver",
+                  "Solver type:\n\t"
+                  "0: Heat\n\t"
+                  "1: P-Laplacian\n\t"
+                  "2: Rvachev scaling");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem type:\n\t"
+                  "0: Point source\n\t"
+                  "1: Circle / sphere level set in 2D / 3D\n\t"
+                  "2: 2D sine-looking level set\n\t"
+                  "3: Gyroid level set in 2D or 3D\n\t"
+                  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&t_param, "-t", "--t-param",
+                  "Diffusion time step (scaled internally scaled by dx*dx).");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0) { args.PrintOptions(cout); }
+
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+
+   // MPI distribution.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   Coefficient *ls_coeff = nullptr;
+   ls_coeff = new FunctionCoefficient(exact_dist_sphere);
+
+   H1_FECollection hfec(order, dim);
+   L2_FECollection lfec(order, dim);
+   ParFiniteElementSpace pfes_s(&pmesh, &hfec), pfes_v(&pmesh, &lfec, dim);
+   ParGridFunction distance_s(&pfes_s), distance_v(&pfes_v);
+
+   distance_s.ProjectCoefficient(*ls_coeff);
+
+
+   ElementMarker marker(pmesh,false, true,2);
+   marker.SetLevelSetFunction(distance_s);
+
+   Array<int> el_marks;
+   marker.MarkElements(el_marks);
+   for(int i=0;i<pmesh.GetNE();i++)
+   {
+       pmesh.SetAttribute(i,el_marks[i]);
+   }
+
+   {
+       GradientGridFunctionCoefficient gco(&distance_s);
+       distance_v.ProjectCoefficient(gco,(int)ElementMarker::SBElementType::CUT);
+   }
+
+
+   ParNonlinearForm* nf=new ParNonlinearForm(&pfes_s);
+   ParNonlinearForm* sf=new ParNonlinearForm(&pfes_s);
+
+   ConstantCoefficient one(1.0);
+   VolShapeIntegrator* itg=new VolShapeIntegrator(one,el_marks);
+   nf->AddDomainIntegrator(itg);
+
+
+   ParGridFunction gyro_s(&pfes_s), gyro_v(&pfes_v);
+   {
+       FunctionCoefficient gfc(Gyroid);
+       gyro_s.ProjectCoefficient(gfc);
+   }
+   GridFunctionCoefficient gfc(&gyro_s);
+   GradientGridFunctionCoefficient gco(&gyro_s);
+
+   SurfShapeIntegrator* its=new SurfShapeIntegrator(gfc, gco, el_marks);
+   sf->AddDomainIntegrator(its);
+
+   Vector lsf(pfes_s.GetTrueVSize());
+   Vector lgr(pfes_s.GetTrueVSize());
+   Vector sgr(pfes_s.GetTrueVSize());
+   distance_s.GetTrueDofs(lsf);
+   double vol=nf->GetEnergy(lsf);
+   double surf=sf->GetEnergy(lsf);
+
+   nf->Mult(lsf,lgr);
+   ParGridFunction grad_s(&pfes_s);
+   grad_s.SetFromTrueDofs(lgr);
+
+   sf->Mult(lsf,sgr);
+
+   if(myrank==0){
+       std::cout<<"Vol="<<vol<<" surf="<<surf<<std::endl;
+   }
+
+
+   //check gradients by FD
+   {
+       mfem::Vector prtv;
+       mfem::Vector tmpv;
+
+       prtv.SetSize(lsf.Size());
+       tmpv.SetSize(lsf.Size());
+
+       prtv.Randomize();
+
+       double nd=mfem::InnerProduct(pmesh.GetComm(),prtv,prtv);
+       double td=mfem::InnerProduct(pmesh.GetComm(),prtv,lgr);
+       double sd=mfem::InnerProduct(pmesh.GetComm(),prtv,sgr);
+
+       td=td/nd;
+       sd=sd/nd;
+       double lsc=1.0;
+       double lqoi;
+       double sqoi;
+
+       for(int l=0; l<8;l++){
+           lsc/=10.0;
+           prtv/=10.0;
+           add(prtv,lsf,tmpv);
+           distance_s.SetFromTrueDofs(tmpv);
+           marker.SetLevelSetFunction(distance_s);
+           marker.MarkElements(el_marks);
+
+           lqoi=nf->GetEnergy(tmpv);
+           sqoi=sf->GetEnergy(tmpv);
+
+           double ld=(lqoi-vol)/lsc;
+           double ls=(sqoi-surf)/lsc;
+           if(myrank==0){
+               std::cout<<" vol="<<vol<<" lvo="<< lqoi<<" dx="<<lsc<<" FD app="<< ld/nd<<" gr="<< td <<" err="<< std::fabs(ld/nd-td) <<std::endl;
+               std::cout<<" sur="<<surf<<" lso="<< sqoi<<" dx="<<lsc<<" FD app="<< ls/nd<<" gr="<< sd<<" err="<< std::fabs(ls/nd-sd) <<std::endl;
+           }
+       }
+
+
+
+   }
+
+   delete sf;
+   delete nf;
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.RegisterField("distance", &distance_s);
+   dacol.RegisterField("grad", &grad_s);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+
+
+   delete ls_coeff;
+
+   return 0;
+}

miniapps/shifted/stress_bracket.cpp

@@ -0,0 +1,221 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+#include "MMA.hpp"
+
+#include "mtop_solvers.hpp"
+#include "mtop_filters.hpp"
+
+class StressBracket2D
+{
+public:
+    StressBracket2D(mfem::ParMesh* mesh_,int order_=1):E(),nu(0.2)
+    {
+        pmesh=mesh_;
+        vorder=order_;
+        esolv=new mfem::ElasticitySolver(pmesh,vorder);
+        esolv->AddMaterial(new mfem::LinIsoElasticityCoefficient(E,nu));
+        esolv->SetNewtonSolver(1e-8,1e-12,1,0);
+        esolv->SetLinearSolver(1e-10,1e-12,400);
+
+        dfes=nullptr;
+        cobj=new mfem::ComplianceObjective();
+
+    }
+
+    ~StressBracket2D()
+    {
+        delete cobj;
+        delete esolv;
+    }
+
+    void SetDesignFES(mfem::ParFiniteElementSpace* fes)
+    {
+        dfes=fes;
+        pdens.SetSpace(dfes);
+    }
+
+    void SetDensity(mfem::Vector& vdens_,
+                    double eta=0.5, double beta=8.0,double pen=3.0){
+        pdens.SetFromTrueDofs(vdens_);
+
+        E.SetDens(&pdens);
+        E.SetProjParam(eta,beta);
+        E.SetEMaxMin(1e-6,1.0);
+        E.SetPenal(pen);
+
+        cobj->SetE(&E);
+        cobj->SetDens(pdens.GetTrueVector());
+        cobj->SetDesignFES(dfes);
+
+        // obtain the response
+        esolv->DelDispBC();
+        esolv->AddDispBC(1,4,0.0);
+        esolv->AddSurfLoad(3,0.00,-1.00,0.0);
+        esolv->FSolve();
+    }
+
+    double Compliance()
+    {
+        return cobj->Eval(esolv->GetDisplacements());
+    }
+
+    double Compliance(mfem::Vector& grad)
+    {
+        cobj->Grad(esolv->GetDisplacements(),grad);
+        return cobj->Eval(esolv->GetDisplacements());
+    }
+
+    mfem::ParGridFunction& GetDisplacements()
+    {
+        return esolv->GetDisplacements();
+    }
+
+
+private:
+    mfem::ParMesh* pmesh;
+
+    mfem::ElasticitySolver* esolv;
+    int vorder;
+
+    mfem::ComplianceObjective* cobj;
+
+    mfem::YoungModulus E;
+    double nu;
+
+    mfem::ParFiniteElementSpace* dfes;
+    mfem::ParGridFunction pdens;
+};
+
+using namespace mfem;
+using namespace std;
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI.
+   int nprocs, myrank;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int solver_type = 0;
+   int rs_levels = 0;
+   int par_ref_levels = 0;
+   int order = 2;
+   int cut_int_order = order;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   double stiff_ratio=1e-6;
+   const char *petscrc_file = "";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&stiff_ratio,"-sr", "--stiff_ratio",
+                  "Stiffness ratio");
+   args.AddOption(&petscrc_file, "-petscopts", "--petscopts",
+                     "PetscOptions file to use.");
+   args.AddOption(&par_ref_levels,
+                     "-rp",
+                     "--refine-parallel",
+                     "Number of times to refine the mesh uniformly in parallel.");
+
+   args.Parse();
+
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0) { args.PrintOptions(cout); }
+
+   mfem::MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.SpaceDimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+   if(myrank==0){
+       std::cout<<"Num elements="<<mesh.GetNE()<<std::endl;
+   }
+
+   // Define a parallel mesh by a partitioning of the serial mesh. Refine
+   // this mesh further in parallel to increase the resolution. Once the
+   // parallel mesh is defined, the serial mesh can be deleted.
+   mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+       for (int l = 0; l < par_ref_levels; l++)
+       {
+          pmesh.UniformRefinement();
+       }
+   }
+
+   if(myrank==0)
+   {
+       std::cout<<"num el="<<pmesh.GetNE()<<std::endl;
+   }
+
+   //allocate the filter
+   mfem::FilterSolver* fsolv=new mfem::FilterSolver(0.02,&pmesh);
+   fsolv->SetSolver(1e-8,1e-12,100,0);
+   fsolv->AddBC(2,0.0);
+   fsolv->AddBC(3,1.0);
+
+   mfem::ParGridFunction pgdens(fsolv->GetFilterFES());
+   mfem::ParGridFunction oddens(fsolv->GetDesignFES());
+   mfem::Vector vdens; vdens.SetSize(fsolv->GetFilterFES()->GetTrueVSize()); vdens=0.0;
+   mfem::Vector vtmpv; vtmpv.SetSize(fsolv->GetDesignFES()->GetTrueVSize()); vtmpv=0.5;
+   oddens.SetFromTrueDofs(vtmpv);
+
+   fsolv->Mult(vtmpv,vdens);
+   pgdens.SetFromTrueDofs(vdens);
+
+   StressBracket2D* brac=new StressBracket2D(&pmesh,1);
+   brac->SetDesignFES(fsolv->GetFilterFES());
+   brac->SetDensity(vdens,0.5,8.0,1.0);
+   brac->Compliance();
+
+
+
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaView", &pmesh);
+   dacol.SetLevelsOfDetail(order);
+   dacol.SetHighOrderOutput(true);
+   dacol.RegisterField("design", &oddens);
+   dacol.RegisterField("flter", &pgdens);
+   dacol.RegisterField("displ", &(brac->GetDisplacements()));
+
+
+
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+
+   delete brac;
+   delete fsolv;
+
+   mfem::MFEMFinalizePetsc();
+   MPI_Finalize();
+   return 0;
+
+}

miniapps/shifted/volp_test.cpp

@@ -0,0 +1,200 @@
+#include <fstream>
+#include <iostream>
+
+#include "mfem.hpp"
+#include "marking.hpp"
+#include "mtop_solvers.hpp"
+#include "mtop_filters.hpp"
+#include "shape_grad.hpp"
+
+using namespace mfem;
+using namespace std;
+
+
+class GyroidCoeff:public Coefficient
+{
+public:
+    GyroidCoeff(double cell_size=1.0){
+        ll=cell_size;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+                const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+        double x = xx[0]*ll;
+        double y = xx[1]*ll;
+        double z = (xx.Size()==3) ? xx[2]*ll : 0.0;
+
+        double r=std::sin(x)*std::cos(y) +
+                std::sin(y)*std::cos(z) +
+               std::sin(z)*std::cos(x) ;
+
+        return r;
+
+        //if(r>0.0){return 1.0;}
+        //return -1.0;
+    }
+
+private:
+    double ll;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int myrank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int solver_type = 0;
+   int rs_levels = 2;
+   int order = 2;
+   const char *device_config = "cpu";
+   bool visualization = true;
+
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myrank == 0) { args.PrintOptions(cout); }
+
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+
+   // MPI distribution.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   int mesh_poly_deg=2;
+   FiniteElementCollection *fec = new H1_FECollection(mesh_poly_deg, dim);
+   ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(&pmesh, fec, dim, Ordering::byNODES);
+
+   pmesh.SetNodalFESpace(pfespace);
+   ParGridFunction x(pfespace);
+   ParGridFunction dx(pfespace);
+   pmesh.SetNodalGridFunction(&x);
+
+   //gradients
+   Vector npos(pfespace->GetTrueVSize());
+   Vector ngrad(pfespace->GetTrueVSize());
+   x.GetTrueDofs(npos);
+
+   GyroidCoeff gyro_co(2*M_PI);
+   ParFiniteElementSpace *cfespace = new ParFiniteElementSpace(&pmesh, fec, 1, Ordering::byNODES);
+   ParGridFunction ggf(cfespace);
+   ggf.ProjectCoefficient(gyro_co);
+
+   ParNonlinearForm* nf=new ParNonlinearForm(pfespace);
+   VolPenalIntegrator* igr=new VolPenalIntegrator(&gyro_co,0.5);
+   nf->AddDomainIntegrator(igr);
+
+   double obj=nf->GetEnergy(npos);
+   std::cout<<"Obj="<<obj<<std::endl;
+
+   nf->Mult(npos,ngrad);
+   dx.SetFromTrueDofs(ngrad);
+
+   // ParaView output.
+   ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
+   dacol.SetLevelsOfDetail(2);
+   dacol.SetHighOrderOutput(true);
+   dacol.RegisterField("coo",&x);
+   dacol.RegisterField("dcoo",&dx);
+   dacol.RegisterField("gyro",&ggf);
+   dacol.SetTime(1.0);
+   dacol.SetCycle(1);
+   dacol.Save();
+
+   //FD check
+   {
+       mfem::Vector prtv;
+       mfem::Vector tmpv;
+
+       prtv.SetSize(npos.Size());
+       tmpv.SetSize(npos.Size());
+
+       prtv.Randomize();
+       double nd=mfem::InnerProduct(pmesh.GetComm(),prtv,prtv);
+       double td=mfem::InnerProduct(pmesh.GetComm(),prtv,ngrad);
+
+       td=td/nd;
+       double lsc=1.0;
+       double lqoi;
+
+       for(int l=0; l<8;l++){
+           lsc/=10.0;
+           prtv/=10.0;
+           add(prtv,npos,tmpv);
+           x.SetFromTrueDofs(tmpv);
+           pmesh.DeleteGeometricFactors();
+           lqoi=nf->GetEnergy(npos);
+           double ld=(lqoi-obj)/lsc;
+           if(myrank==0){
+               std::cout<<" obj="<<obj<<" lvo="<< lqoi<<" dx="<<lsc<<" FD app="<< ld/nd<<" gr="<< td <<" err="<< std::fabs(ld/nd-td) <<std::endl;
+           }
+       }
+   }
+
+   {
+       mfem::Vector tmpv(npos);
+       tmpv.SetSize(npos.Size());
+       double lqoi;
+       for(int i=0;i<800;i++){
+          tmpv.Add(-0.0005,ngrad);
+          x.SetFromTrueDofs(tmpv);
+          pmesh.DeleteGeometricFactors();
+          lqoi=nf->GetEnergy(npos);
+          nf->Mult(npos,ngrad);
+          std::cout<<"lqoi="<<lqoi<<" nr="<<ngrad.Norml2()<<std::endl;
+       }
+   }
+
+   ggf.ProjectCoefficient(gyro_co);
+
+
+   dacol.SetTime(2.0);
+   dacol.SetCycle(2);
+   dacol.Save();
+
+
+
+
+   delete nf;
+   delete cfespace;
+   delete pfespace;
+   delete fec;
+
+   return 0;
+
+}