v0.9.1
basic_helmholtz.cpp

Using Basic interface calculate the divergence of base functions, and integral of flux on the boundary. Since the h-div space is used, volume integral and boundary integral should give the same result.

/**
* \file basic_helmholtz.cpp
* \example basic_helmholtz.cpp
*
* Using Basic interface calculate the divergence of base functions, and
* integral of flux on the boundary. Since the h-div space is used, volume
* integral and boundary integral should give the same result.
*/
/* This file is part of MoFEM.
* MoFEM 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 3 of the License, or (at your
* option) any later version.
*
* MoFEM 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. 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 MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
#include <BaseOps.hpp>
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode setUP();
MoFEMErrorCode kspSolve();
MoFEMErrorCode postProcess();
MatrixDouble invJac;
};
CHKERR setUP();
CHKERR bC();
CHKERR OPs();
CHKERR kspSolve();
CHKERR postProcess();
}
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// Add field
CHKERR simple->addDomainField("U_REAL", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
1);
CHKERR simple->addDomainField("U_IMAG", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
1);
CHKERR simple->addBoundaryField("U_REAL", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
1);
CHKERR simple->addBoundaryField("U_IMAG", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
1);
constexpr int order = 10;
CHKERR simple->setFieldOrder("U_REAL", order);
CHKERR simple->setFieldOrder("U_IMAG", order);
CHKERR simple->setUp();
}
//! [Set up problem]
//! [Applying essential BC]
}
//! [Applying essential BC]
//! [Push operators to pipeline]
Basic *basic = mField.getInterface<Basic>();
auto vol_source_function = [](const double x, const double y,
const double z) {
const double xc = 0;
const double yc = -1.25;
const double xs = x - xc;
const double ys = y - yc;
constexpr double eps = 60;
return exp(-pow(eps * sqrt(xs * xs + ys * ys), 2));
};
constexpr double k = 90;
auto beta = [](const double, const double, const double) { return -1; };
auto k2 = [k](const double, const double, const double) { return pow(k, 2); };
auto kp = [k](const double, const double, const double) { return k; };
auto km = [k](const double, const double, const double) { return -k; };
auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
auto set_domain = [&]() {
basic->getOpDomainLhsPipeline().push_back(
basic->getOpDomainLhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
basic->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("U_REAL", "U_REAL", beta));
basic->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("U_IMAG", "U_IMAG", beta));
basic->getOpDomainLhsPipeline().push_back(
new OpDomainMass("U_REAL", "U_REAL", k2));
basic->getOpDomainLhsPipeline().push_back(
new OpDomainMass("U_IMAG", "U_IMAG", k2));
basic->getOpDomainRhsPipeline().push_back(
new OpDomainSource("U_REAL", vol_source_function));
CHKERR basic->setDomainRhsIntegrationRule(integration_rule);
CHKERR basic->setDomainLhsIntegrationRule(integration_rule);
};
auto set_boundary = [&]() {
basic->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass("U_IMAG", "U_REAL", kp));
basic->getOpBoundaryLhsPipeline().push_back(
new OpBoundaryMass("U_REAL", "U_IMAG", km));
CHKERR basic->setBoundaryLhsIntegrationRule(integration_rule);
};
CHKERR set_domain();
CHKERR set_boundary();
}
//! [Push operators to pipeline]
//! [Solve]
Simple *simple = mField.getInterface<Simple>();
Basic *basic = mField.getInterface<Basic>();
auto solver = basic->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto dm = simple->getDM();
auto D = smartCreateDMDVector(dm);
CHKERR KSPSolve(solver, F, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
}
//! [Solve]
//! [Postprocess results]
Basic *basic = mField.getInterface<Basic>();
basic->getDomainLhsFE().reset();
basic->getDomainRhsFE().reset();
basic->getBoundaryLhsFE().reset();
basic->getBoundaryRhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->addFieldValuesPostProc("U_REAL");
post_proc_fe->addFieldValuesPostProc("U_IMAG");
basic->getDomainRhsFE() = post_proc_fe;
CHKERR post_proc_fe->writeFile("out_helmholtz.h5m");
}
//! [Postprocess results]
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Load mesh]
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Load mesh]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}