13 static char help[] =
"...\n\n";
35 return 200 * sin(x * 10.) * cos(y * 10.);
41 return sin(x * 10.) * cos(y * 10.);
63 : domainField(
"U"), boundaryField(
"L"), mField(m_field) {}
97 auto get_ents_on_mesh = [&]() {
98 Range boundary_entities;
100 std::string entity_name = it->getName();
101 if (entity_name.compare(0, 18,
"BOUNDARY_CONDITION") == 0) {
103 boundary_entities,
true);
107 Range boundary_vertices;
109 boundary_vertices,
true);
110 boundary_entities.merge(boundary_vertices);
115 return boundary_entities;
118 auto mark_boundary_dofs = [&](
Range &&skin_edges) {
120 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
122 skin_edges, *marker_ptr);
138 auto det_ptr = boost::make_shared<VectorDouble>();
139 auto jac_ptr = boost::make_shared<MatrixDouble>();
140 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
144 pipeline_mng->getOpDomainLhsPipeline().push_back(
146 pipeline_mng->getOpDomainLhsPipeline().push_back(
148 pipeline_mng->getOpDomainLhsPipeline().push_back(
150 pipeline_mng->getOpDomainLhsPipeline().push_back(
152 pipeline_mng->getOpDomainLhsPipeline().push_back(
158 pipeline_mng->getOpDomainRhsPipeline().push_back(
164 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
166 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
172 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
184 auto domain_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * (p - 1); };
185 auto domain_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * (p - 1); };
186 CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
187 CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
189 auto boundary_rule_lhs = [](
int,
int,
int p) ->
int {
return 2 * p; };
190 auto boundary_rule_rhs = [](
int,
int,
int p) ->
int {
return 2 * p; };
191 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
192 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
202 auto ksp_solver = pipeline_mng->
createKSP();
203 CHKERR KSPSetFromOptions(ksp_solver);
213 CHKERR KSPGetPC(ksp_solver, &pc);
214 PetscBool is_pcfs = PETSC_FALSE;
215 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
220 if (is_pcfs == PETSC_TRUE) {
221 IS is_domain, is_boundary;
222 cerr <<
"Running FIELDSPLIT..." << endl;
230 CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
232 CHKERR ISDestroy(&is_boundary);
236 CHKERR KSPSetUp(ksp_solver);
242 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
243 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
254 pipeline_mng->getBoundaryLhsFE().reset();
256 auto d_ptr = boost::make_shared<VectorDouble>();
257 auto l_ptr = boost::make_shared<VectorDouble>();
261 auto post_proc_domain_fe = boost::make_shared<PostProcFaceEle>(
mField);
262 post_proc_domain_fe->getOpPtrVector().push_back(
264 post_proc_domain_fe->getOpPtrVector().push_back(
265 new OpPPMap(post_proc_domain_fe->getPostProcMesh(),
266 post_proc_domain_fe->getMapGaussPts(), {{domainField, d_ptr}},
268 pipeline_mng->getDomainRhsFE() = post_proc_domain_fe;
270 auto post_proc_boundary_fe = boost::make_shared<PostProcEdgeEle>(mField);
271 post_proc_boundary_fe->getOpPtrVector().push_back(
273 post_proc_boundary_fe->getOpPtrVector().push_back(
274 new OpPPMap(post_proc_boundary_fe->getPostProcMesh(),
275 post_proc_boundary_fe->getMapGaussPts(),
276 {{boundaryField, l_ptr}}, {}, {}, {}));
277 pipeline_mng->getBoundaryRhsFE() = post_proc_boundary_fe;
279 CHKERR pipeline_mng->loopFiniteElements();
280 CHKERR post_proc_domain_fe->writeFile(
"out_result_domain.h5m");
281 CHKERR post_proc_boundary_fe->writeFile(
"out_result_boundary.h5m");
300 int main(
int argc,
char *argv[]) {
303 const char param_file[] =
"param_file.petsc";
309 DMType dm_name =
"DMMOFEM";