#ifndef __POISSON2DLAGRANGEMULTIPLIER_HPP__
#define __POISSON2DLAGRANGEMULTIPLIER_HPP__
#include <stdlib.h>
typedef boost::function<
double(
const double,
const double,
const double)>
public:
OpDomainLhsK(std::string row_field_name, std::string col_field_name,
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
boundaryMarker(boundary_marker) {
sYmm = true;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locLhs.resize(nb_row_dofs, nb_col_dofs, false);
locLhs.clear();
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
locLhs(rr, cc) += t_row_diff_base(
i) * t_col_diff_base(
i) *
a;
++t_col_diff_base;
}
++t_row_diff_base;
}
++t_w;
}
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transLocLhs) = trans(locLhs);
ADD_VALUES);
}
}
}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
public:
boost::shared_ptr<std::vector<unsigned char>> boundary_marker = nullptr)
sourceTermFunc(source_term_function), boundaryMarker(boundary_marker) {}
if (nb_dofs) {
locRhs.resize(nb_dofs, false);
locRhs.clear();
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
double body_source =
sourceTermFunc(t_coords(0), t_coords(1), t_coords(2));
for (int rr = 0; rr != nb_dofs; rr++) {
locRhs[rr] += t_base * body_source *
a;
++t_base;
}
++t_w;
++t_coords;
}
}
}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
};
public:
OpBoundaryLhsC(std::string row_field_name, std::string col_field_name)
sYmm = false;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
locBoundaryLhs.resize(nb_row_dofs, nb_col_dofs, false);
locBoundaryLhs.clear();
const double edge = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; cc++) {
locBoundaryLhs(rr, cc) += t_row_base * t_col_base *
a;
++t_col_base;
}
++t_row_base;
}
++t_w;
}
ADD_VALUES);
}
}
private:
};
public:
boundaryFunc(boundary_function) {}
if (nb_dofs) {
locBoundaryRhs.resize(nb_dofs, false);
locBoundaryRhs.clear();
const double edge = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * edge;
double boundary_term =
boundaryFunc(t_coords(0), t_coords(1), t_coords(2));
for (int rr = 0; rr != nb_dofs; rr++) {
locBoundaryRhs[rr] += t_base * boundary_term *
a;
++t_base;
}
++t_w;
++t_coords;
}
ADD_VALUES);
}
}
private:
};
};
#endif //__POISSON2DLAGRANGEMULTIPLIER_HPP__
#include <stdlib.h>
static char help[] =
"...\n\n";
public:
private:
static double sourceTermFunction(const double x, const double y,
const double z) {
return 200 * sin(x * 10.) * cos(y * 10.);
}
static double boundaryFunction(const double x, const double y,
const double z) {
return sin(x * 10.) * cos(y * 10.);
}
std::string boundaryField;
int oRder;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
Range boundaryEntitiesForFieldsplit;
};
:
domainField(
"U"), boundaryField(
"L"), mField(m_field) {}
}
}
auto get_ents_on_mesh = [&]() {
std::string entity_name = it->getName();
if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
boundary_entities, true);
}
}
boundary_vertices, true);
boundary_entities.merge(boundary_vertices);
return boundary_entities;
};
auto mark_boundary_dofs = [&](
Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
{
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
}
{
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
{
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
}
{
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
}
}
auto domain_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * (p - 1); };
auto domain_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * (p - 1); };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
auto boundary_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * p; };
auto boundary_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * p; };
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
}
auto ksp_solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(ksp_solver);
if (1) {
PC pc;
CHKERR KSPGetPC(ksp_solver, &pc);
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
IS is_domain, is_boundary;
cerr << "Running FIELDSPLIT..." << endl;
CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
CHKERR ISDestroy(&is_boundary);
}
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto d_ptr = boost::make_shared<VectorDouble>();
auto l_ptr = boost::make_shared<VectorDouble>();
auto post_proc_domain_fe = boost::make_shared<PostProcFaceEle>(
mField);
post_proc_domain_fe->getOpPtrVector().push_back(
post_proc_domain_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_domain_fe->getPostProcMesh(),
post_proc_domain_fe->getMapGaussPts(), {{domainField, d_ptr}},
{}, {}, {}));
pipeline_mng->getDomainRhsFE() = post_proc_domain_fe;
auto post_proc_boundary_fe = boost::make_shared<PostProcEdgeEle>(mField);
post_proc_boundary_fe->getOpPtrVector().push_back(
post_proc_boundary_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_boundary_fe->getPostProcMesh(),
post_proc_boundary_fe->getMapGaussPts(),
{{boundaryField, l_ptr}}, {}, {}, {}));
pipeline_mng->getBoundaryRhsFE() = post_proc_boundary_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_domain_fe->writeFile(
"out_result_domain.h5m");
CHKERR post_proc_boundary_fe->writeFile(
"out_result_boundary.h5m");
}
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
try {
DMType dm_name = "DMMOFEM";
CHKERR poisson_problem.runProgram();
}
return 0;
}