diffusion-3d.edp

/*
   Author(s): Pierre Jolivet <pierre@joliv.et>
        Date: 2019-02-18

   Copyright (C) 2019-     Centre National de la Recherche Scientifique

   This script is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

   If you use this script, you are kindly asked to cite the following article:

   "A Multilevel Schwarz Preconditioner Based on a Hierarchy of Robust Coarse Spaces",
   H. Al Daas, L. Grigori, P. Jolivet, and P.-H. Tournier (2019).
 */

IFMACRO(petsc)
load "PETSc"
bool withPETSc = true;
ENDIFMACRO
IFMACRO(!petsc)
load "hpddm"
bool withPETSc = false;
ENDIFMACRO

macro dimension()3// EOM
macro def(i)i// EOM
macro init(i)i// EOM
include "coefficients.idp"
func Pk = P2;                     // finite element space
int[int] LL = [2, 3, 2, 1, 2, 2];
include "macro_ddm.idp"           // additional DDM functions

meshN Th = cube(getARGV("-global", 10), getARGV("-global", 10), getARGV("-global", 10), label = LL);
fespace Wh(Th, Pk);               // local finite element space
int[int][int] intersection;       // local-to-neighbors renumbering
real[int] D;                      // partition of unity
{
    int s = getARGV("-split", 1); // refinement factor
    build(Th, s, intersection, D, Pk, mpiCommWorld)
}

real[int] rhs(Wh.ndof);           // local right-hand side
matrix<real> Loc, Neumann;        // local operator
Wh<real> def(Rb)[0];              // rigid body modes
{                                 // local weak form
    fespace Ph(Th, P0);
    Ph kappa = skyscraper(x, y);
    real f = -10000.0;
    varf vPb(def(u), def(v)) = intN(Th, qforder = 3)(kappa * grad(u)' * grad(v)) + intN(Th)(f * v) + on(1, u = 0.0);
    Loc = vPb(Wh, Wh, sym = !withPETSc, tgv = -2);
    rhs = vPb(0, Wh, tgv = -1);
}

macro prefix()"diffusion_3d"// EOM
include "common.idp"