v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
SetUpSchurImpl Struct Reference
Inheritance diagram for SetUpSchurImpl:
[legend]
Collaboration diagram for SetUpSchurImpl:
[legend]

Public Member Functions

 SetUpSchurImpl (MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (KSP solver)
 
MoFEMErrorCode preProc ()
 
MoFEMErrorCode postProc ()
 
 SetUpSchurImpl (MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_is, SmartPetscObj< AO > ao_up)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (TS solver)
 
MoFEMErrorCode preProc ()
 
MoFEMErrorCode postProc ()
 
 SetUpSchurImpl (MoFEM::Interface &m_field)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (SmartPetscObj< TS > solver)
 
 SetUpSchurImpl (MoFEM::Interface &m_field)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (SmartPetscObj< TS > solver)
 
 SetUpSchurImpl (MoFEM::Interface &m_field)
 
virtual ~SetUpSchurImpl ()
 
MoFEMErrorCode setUp (SmartPetscObj< KSP > solver)
 
virtual MoFEMErrorCode setUp (TS solver)=0
 
virtual MoFEMErrorCode setUp (SmartPetscObj< TS > solver)=0
 
virtual MoFEMErrorCode setUp (SmartPetscObj< TS > solver)=0
 
virtual MoFEMErrorCode setUp (SmartPetscObj< KSP > solver)=0
 

Private Member Functions

MoFEMErrorCode setEntities ()
 
MoFEMErrorCode setOperator ()
 
MoFEMErrorCode setPC (PC pc)
 
SmartPetscObj< DM > createSubDM ()
 
MoFEMErrorCode setOperator ()
 
MoFEMErrorCode setPC (PC pc)
 
SmartPetscObj< DM > createSubDM ()
 
MoFEMErrorCode setEntities ()
 
MoFEMErrorCode setUpSubDM ()
 
MoFEMErrorCode setOperator ()
 
MoFEMErrorCode setPC (PC pc)
 

Private Attributes

SmartPetscObj< Mat > M
 
SmartPetscObj< Mat > P
 
SmartPetscObj< Mat > S
 
SmartPetscObj< AO > aoUp
 
SmartPetscObj< IS > isA00
 
SmartPetscObj< IS > isShift
 
MoFEM::InterfacemField
 
EshelbianCore * epCorePtr
 
SmartPetscObj< DM > subDM
 field split sub dm More...
 
SmartPetscObj< IS > fieldSplitIS
 IS for split Schur block. More...
 
SmartPetscObj< Mat > A
 
Range volEnts
 
Range subEnts
 

Additional Inherited Members

- Static Public Member Functions inherited from SetUpSchur
static boost::shared_ptr< SetUpSchurcreateSetUpSchur (MoFEM::Interface &m_field, SmartPetscObj< DM > sub_dm, SmartPetscObj< IS > field_split_it, SmartPetscObj< AO > ao_map)
 Create data structure for handling Schur complement. More...
 
static boost::shared_ptr< SetUpSchurcreateSetUpSchur (MoFEM::Interface &m_field)
 
static boost::shared_ptr< SetUpSchurcreateSetUpSchur (MoFEM::Interface &m_field)
 
static boost::shared_ptr< SetUpSchurcreateSetUpSchur (MoFEM::Interface &m_field)
 
- Protected Member Functions inherited from SetUpSchur
 SetUpSchur ()=default
 
 SetUpSchur ()=default
 
 SetUpSchur ()=default
 
 SetUpSchur ()=default
 

Detailed Description

Examples
plastic.cpp.

Definition at line 9 of file SetUpSchurImpl.cpp.

Constructor & Destructor Documentation

◆ SetUpSchurImpl() [1/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field,
SmartPetscObj< Mat >  m,
SmartPetscObj< Mat >  p,
EshelbianCore *  ep_core_ptr 
)
inline

Definition at line 11 of file SetUpSchurImpl.cpp.

13 : SetUpSchur(), mField(m_field), M(m), P(p), epCorePtr(ep_core_ptr) {
14 if (S) {
17 "Is expected that schur matrix is not allocated. This is "
18 "possible only is PC is set up twice");
19 }
20 }
static Index< 'p', 3 > p
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FTensor::Index< 'm', SPACE_DIM > m
SetUpSchur()=default
EshelbianCore * epCorePtr
SmartPetscObj< Mat > S
SmartPetscObj< Mat > P
SmartPetscObj< Mat > M
MoFEM::Interface & mField

◆ ~SetUpSchurImpl() [1/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual
Examples
plastic.cpp.

Definition at line 21 of file SetUpSchurImpl.cpp.

21{ S.reset(); }

◆ SetUpSchurImpl() [2/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field,
SmartPetscObj< DM >  sub_dm,
SmartPetscObj< IS >  field_split_is,
SmartPetscObj< AO >  ao_up 
)
inline

Definition at line 1273 of file plastic.cpp.

1275 : SetUpSchur(), mField(m_field), subDM(sub_dm),
1276 fieldSplitIS(field_split_is), aoUp(ao_up) {
1277 if (S) {
1280 "Is expected that schur matrix is not allocated. This is "
1281 "possible only is PC is set up twice");
1282 }
1283 }
SmartPetscObj< DM > subDM
field split sub dm
Definition: plastic.cpp:1304
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
Definition: plastic.cpp:1305
SmartPetscObj< AO > aoUp

◆ ~SetUpSchurImpl() [2/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 1284 of file plastic.cpp.

1284 {
1285#ifdef ADD_CONTACT
1286 A.reset();
1287 P.reset();
1288#endif // ADD_CONTACT
1289 S.reset();
1290 }
SmartPetscObj< Mat > A
Definition: contact.cpp:983

◆ SetUpSchurImpl() [3/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field)
inline

Definition at line 966 of file contact.cpp.

966: SetUpSchur(), mField(m_field) {}

◆ ~SetUpSchurImpl() [3/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 968 of file contact.cpp.

968 {
969 A.reset();
970 P.reset();
971 S.reset();
972 }

◆ SetUpSchurImpl() [4/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field)
inline

Definition at line 935 of file incompressible_elasticity.cpp.

935: SetUpSchur(), mField(m_field) {}

◆ ~SetUpSchurImpl() [4/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 936 of file incompressible_elasticity.cpp.

936{ S.reset(); }

◆ SetUpSchurImpl() [5/5]

SetUpSchurImpl::SetUpSchurImpl ( MoFEM::Interface m_field)
inline

Definition at line 883 of file elastic.cpp.

883 : SetUpSchur(), mField(m_field) {
884 if (S) {
887 "Is expected that schur matrix is not allocated. This is "
888 "possible only is PC is set up twice");
889 }
890 }

◆ ~SetUpSchurImpl() [5/5]

virtual SetUpSchurImpl::~SetUpSchurImpl ( )
inlinevirtual

Definition at line 891 of file elastic.cpp.

891{ S.reset(); }

Member Function Documentation

◆ createSubDM() [1/2]

SmartPetscObj< DM > SetUpSchurImpl::createSubDM ( )
private

Definition at line 1056 of file contact.cpp.

1056 {
1057 auto simple = mField.getInterface<Simple>();
1058 auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
1059 auto set_up = [&]() {
1061 CHKERR DMMoFEMCreateSubDM(sub_dm, simple->getDM(), "SUB");
1062 CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
1063 CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
1064 CHKERR DMMoFEMAddSubFieldRow(sub_dm, "U");
1065 CHKERR DMSetUp(sub_dm);
1067 };
1068 CHK_THROW_MESSAGE(set_up(), "sey up dm");
1069 return sub_dm;
1070}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:219
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:442
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
virtual MPI_Comm & get_comm() const =0
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ createSubDM() [2/2]

SmartPetscObj< DM > SetUpSchurImpl::createSubDM ( )
private

◆ postProc() [1/2]

MoFEMErrorCode SetUpSchurImpl::postProc ( )
Examples
plastic.cpp.

Definition at line 204 of file SetUpSchurImpl.cpp.

204 {
206 if (S) {
207 MOFEM_LOG("EP", Sev::verbose) << "Assemble Schur";
208 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
209 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
210 }
211 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
212 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
213 CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
214 CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
215 CHKERR MatAXPY(P, 1, M, SAME_NONZERO_PATTERN);
216
218}
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308

◆ postProc() [2/2]

MoFEMErrorCode SetUpSchurImpl::postProc ( )

◆ preProc() [1/2]

MoFEMErrorCode SetUpSchurImpl::preProc ( )
Examples
plastic.cpp.

Definition at line 193 of file SetUpSchurImpl.cpp.

193 {
195 MOFEM_LOG("EP", Sev::verbose) << "Zero Schur";
196 if (SetUpSchurImpl::S) {
197 CHKERR MatZeroEntries(S);
198 }
199 CHKERR MatZeroEntries(M);
200 CHKERR MatZeroEntries(P);
202}

◆ preProc() [2/2]

MoFEMErrorCode SetUpSchurImpl::preProc ( )

◆ setEntities() [1/2]

MoFEMErrorCode SetUpSchurImpl::setEntities ( )
private

Definition at line 942 of file elastic.cpp.

942 {
945 CHKERR mField.get_moab().get_entities_by_dimension(simple->getMeshset(),
947 CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
948 subEnts);
949 subEnts = subtract(subEnts, volEnts);
951};
constexpr int SPACE_DIM
[Define dimension]
Definition: elastic.cpp:18
virtual moab::Interface & get_moab()=0

◆ setEntities() [2/2]

MoFEMErrorCode SetUpSchurImpl::setEntities ( )
private

◆ setOperator() [1/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )
private

Definition at line 1072 of file contact.cpp.

1072 {
1074
1075 double eps_stab = 1e-4;
1076 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1077 PETSC_NULL);
1078
1079 using B =
1081 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1082
1083 auto pip = mField.getInterface<PipelineManager>();
1084 // Boundary
1085 auto dm_is = getDMSubData(subDM)->getSmartRowIs();
1086 auto ao_up = createAOMappingIS(dm_is, PETSC_NULL);
1087
1088 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1089 pip->getOpBoundaryLhsPipeline().push_back(
1090 new OpMassStab("SIGMA", "SIGMA",
1091 [eps_stab](double, double, double) { return eps_stab; }));
1092 pip->getOpBoundaryLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
1093 {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
1094
1095 // Domain
1096 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1097 pip->getOpDomainLhsPipeline().push_back(new OpSchurAssembleEnd<SCHUR_DGESV>(
1098 {"SIGMA"}, {nullptr}, {ao_up}, {S}, {false}, false));
1099
1100 auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1101 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1102
1103 pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
1105 CHKERR MatZeroEntries(A);
1106 CHKERR MatZeroEntries(P);
1107 CHKERR MatZeroEntries(S);
1108 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Begin";
1110 };
1111
1112 post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
1113 post_proc_schur_lhs_ptr]() {
1115 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble End";
1116
1117 *post_proc_schur_lhs_ptr->matAssembleSwitch = false;
1118
1119 auto print_mat_norm = [this](auto a, std::string prefix) {
1121 double nrm;
1122 CHKERR MatNorm(a, NORM_FROBENIUS, &nrm);
1123 MOFEM_LOG("CONTACT", Sev::noisy) << prefix << " norm = " << nrm;
1125 };
1126
1127 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1128 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1130 mField, post_proc_schur_lhs_ptr, 1, A)();
1131
1132 CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1133 CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1134 CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
1135
1136 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1137 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1138
1140 mField, post_proc_schur_lhs_ptr, 1, S, ao_up)();
1141
1142#ifndef NDEBUG
1143 CHKERR print_mat_norm(A, "A");
1144 CHKERR print_mat_norm(P, "P");
1145 CHKERR print_mat_norm(S, "S");
1146#endif // NDEBUG
1147
1148 MOFEM_LOG("CONTACT", Sev::verbose) << "Lhs Assemble Finish";
1150 };
1151
1152 auto simple = mField.getInterface<Simple>();
1153 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1154 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
1155 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1156
1158}
constexpr double a
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1061
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:33
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
PipelineManager interface.

◆ setOperator() [2/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )
private

◆ setOperator() [3/3]

MoFEMErrorCode SetUpSchurImpl::setOperator ( )
private

◆ setPC() [1/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc)
private

Definition at line 1160 of file contact.cpp.

1160 {
1162 auto simple = mField.getInterface<Simple>();
1164 mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1165 simple->getProblemName(), ROW, "SIGMA", 0, SPACE_DIM, is);
1166 CHKERR PCFieldSplitSetIS(pc, NULL, is);
1167 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1169}
@ ROW
Definition: definitions.h:123
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
intrusive_ptr for managing petsc objects

◆ setPC() [2/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc)
private

◆ setPC() [3/3]

MoFEMErrorCode SetUpSchurImpl::setPC ( PC  pc)
private

◆ setUp() [1/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( KSP  solver)
Examples
plastic.cpp.

Definition at line 39 of file SetUpSchurImpl.cpp.

39 {
41
42 auto get_ents_by_dim = [&](int dim) {
43 Range ents;
44 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents);
45 return ents;
46 };
47
48 std::vector<std::string> schur_field_list{
49
50 epCorePtr->stretchTensor, epCorePtr->bubbleField,
51
52 epCorePtr->rotAxis, epCorePtr->spatialL2Disp,
53
54 epCorePtr->piolaStress
55
56 };
57 std::vector<std::string> field_list{epCorePtr->piolaStress};
58
59 auto create_schur_dm = [&](SmartPetscObj<DM> &dm_sub) {
61
62 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
63 CHKERR DMMoFEMCreateSubDM(dm_sub, epCorePtr->dmElastic, "SUB_SCHUR");
64 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
65 CHKERR DMMoFEMSetIsPartitioned(dm_sub, PETSC_TRUE);
66 CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->elementVolumeName);
67 CHKERR DMMoFEMAddElement(dm_sub, epCorePtr->naturalBcElement);
68
69 auto faces_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM - 1));
70 std::vector<boost::shared_ptr<Range>> dm_range_list{faces_ptr};
71
72 if (field_list.size() != dm_range_list.size()) {
73 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
74 "Both ranges should have the same size");
75 }
76
77 int r_idx = 0;
78 for (auto f : field_list) {
79 MOFEM_LOG("EP", Sev::inform) << "Add schur field: " << f;
80 CHKERR DMMoFEMAddSubFieldRow(dm_sub, f.c_str(), dm_range_list[r_idx]);
81 CHKERR DMMoFEMAddSubFieldCol(dm_sub, f.c_str(), dm_range_list[r_idx]);
82 ++r_idx;
83 }
84 CHKERR DMSetUp(dm_sub);
86 };
87
88 auto create_a00_is = [&](SmartPetscObj<IS> &is_a00) {
90 auto vols = get_ents_by_dim(SPACE_DIM);
91 std::vector<SmartPetscObj<IS>> is_vec;
92 std::vector<Range *> range_list_ptr(schur_field_list.size(), &vols);
93 int r_idx = 0;
94 for (auto f : schur_field_list) {
95 SmartPetscObj<IS> is;
96 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
97 "ELASTIC_PROBLEM", ROW, f, 0, MAX_DOFS_ON_ENTITY, is,
98 range_list_ptr[r_idx]);
99 is_vec.push_back(is);
100 ++r_idx;
101 }
102 if (is_vec.size()) {
103 is_a00 = is_vec[0];
104 for (auto i = 1; i < is_vec.size(); ++i) {
105 IS is_union_raw;
106 CHKERR ISExpand(is_a00, is_vec[i], &is_union_raw);
107 is_a00 = SmartPetscObj<IS>(is_union_raw);
108 }
109 CHKERR ISSort(is_a00);
110 }
112 };
113
114 PC pc;
115 CHKERR KSPSetFromOptions(solver);
116 CHKERR KSPGetPC(solver, &pc);
117 PetscBool is_pcfs = PETSC_FALSE;
118 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
119 if (is_pcfs) {
120
121 SmartPetscObj<DM> schur_dm, a00_dm;
122 CHKERR create_schur_dm(schur_dm);
123 CHKERR create_a00_is(isA00);
124
125 if (S) {
128 "Is expected that schur matrix is not allocated. This is "
129 "possible only is PC is set up twice");
130 }
131
132 S = createDMMatrix(schur_dm);
133
134 auto set_ops = [&]() {
136 auto range_ptr = boost::make_shared<Range>(get_ents_by_dim(SPACE_DIM));
137 std::vector<boost::shared_ptr<Range>> ranges_list(schur_field_list.size(),
138 range_ptr);
139 std::vector<SmartPetscObj<AO>> ao_list(schur_field_list.size(),
140 SmartPetscObj<AO>());
141 std::vector<SmartPetscObj<Mat>> mat_list(schur_field_list.size(),
142 SmartPetscObj<Mat>());
143 std::vector<bool> symm_list(schur_field_list.size(), false);
144
145 auto dm_is = getDMSubData(schur_dm)->getSmartRowIs();
146 aoUp = createAOMappingIS(dm_is, PETSC_NULL);
147 ao_list.back() = aoUp;
148 mat_list.back() = S;
149
150 const bool symmetric_system =
151 EshelbianCore::gradApperoximator <= MODERATE_ROT &&
152 EshelbianCore::rotSelector == SMALL_ROT;
153 epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
154 new OpSchurAssembleBegin());
155 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
156 new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
157 ao_list, mat_list, symm_list,
158 symmetric_system));
159 epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
160 new OpSchurAssembleBegin());
161 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
162 new OpSchurAssembleEnd<SCHUR_DGESV>(schur_field_list, ranges_list,
163 ao_list, mat_list, symm_list,
164 symmetric_system));
165
167 };
168
169 auto set_pc = [&]() {
171 CHKERR PCFieldSplitSetIS(pc, NULL, isA00);
172 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
174 };
175
176 CHKERR set_ops();
177 CHKERR set_pc();
178
179 } else {
180 epCorePtr->elasticFeLhs->getOpPtrVector().push_front(
181 new OpSchurAssembleBegin());
182 epCorePtr->elasticFeLhs->getOpPtrVector().push_back(
183 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
184 epCorePtr->elasticBcLhs->getOpPtrVector().push_front(
185 new OpSchurAssembleBegin());
186 epCorePtr->elasticBcLhs->getOpPtrVector().push_back(
187 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
188 }
189
191}
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
const int dim
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:275
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:988
FTensor::Index< 'i', SPACE_DIM > i
SmartPetscObj< IS > isA00

◆ setUp() [2/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< KSP >  solver)
virtual

Implements SetUpSchur.

Definition at line 910 of file elastic.cpp.

910 {
912 auto pip = mField.getInterface<PipelineManager>();
913 PC pc;
914 CHKERR KSPGetPC(solver, &pc);
915 PetscBool is_pcfs = PETSC_FALSE;
916 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
917 if (is_pcfs) {
918 if (S) {
921 "Is expected that schur matrix is not allocated. This is "
922 "possible only is PC is set up twice");
923 }
927 CHKERR MatSetBlockSize(S, SPACE_DIM);
928 CHKERR MatSetOption(S, MAT_SYMMETRIC, PETSC_TRUE);
930 CHKERR setPC(pc);
931 } else {
932 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
933 pip->getOpBoundaryLhsPipeline().push_back(
934 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
935 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
936 pip->getOpDomainLhsPipeline().push_back(
937 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
938 }
940}
MoFEMErrorCode setEntities()
Definition: elastic.cpp:942
MoFEMErrorCode setPC(PC pc)
Definition: contact.cpp:1160
MoFEMErrorCode setUpSubDM()
Definition: elastic.cpp:953
MoFEMErrorCode setOperator()
Definition: contact.cpp:1072

◆ setUp() [3/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< TS >  solver)
virtual

Implements SetUpSchur.

Definition at line 992 of file contact.cpp.

992 {
995 auto pip = mField.getInterface<PipelineManager>();
996 auto dm = simple->getDM();
997
998 SNES snes;
999 CHKERR TSGetSNES(solver, &snes);
1000 KSP ksp;
1001 CHKERR SNESGetKSP(snes, &ksp);
1002 CHKERR KSPSetFromOptions(ksp);
1003
1004 PC pc;
1005 CHKERR KSPGetPC(ksp, &pc);
1006
1007 PetscBool is_pcfs = PETSC_FALSE;
1008 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1009 if (is_pcfs) {
1010
1011 MOFEM_LOG("CONTACT", Sev::inform) << "Setup Schur pc";
1012
1013 if (A || P || S) {
1016 "It is expected that Schur matrix is not allocated. This is "
1017 "possible only if PC is set up twice");
1018 }
1019
1020 A = createDMMatrix(dm);
1021 P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
1022 subDM = createSubDM();
1024 CHKERR MatSetBlockSize(S, SPACE_DIM);
1025
1026 auto ts_ctx_ptr = getDMTsCtx(dm);
1027 CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
1028
1030 CHKERR setPC(pc);
1031
1032 } else {
1033 MOFEM_LOG("CONTACT", Sev::inform) << "No Schur pc";
1034
1035 pip->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1036 pip->getOpBoundaryLhsPipeline().push_back(
1037 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1038 pip->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1039 pip->getOpDomainLhsPipeline().push_back(
1040 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1041
1042 auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
1043 post_proc_schur_lhs_ptr->postProcessHook = [this,
1044 post_proc_schur_lhs_ptr]() {
1047 mField, post_proc_schur_lhs_ptr, 1.)();
1049 };
1050 auto ts_ctx_ptr = getDMTsCtx(dm);
1051 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
1052 }
1054}
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
SmartPetscObj< DM > createSubDM()
Definition: contact.cpp:1056

◆ setUp() [4/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( SmartPetscObj< TS >  solver)
virtual

Implements SetUpSchur.

◆ setUp() [5/5]

MoFEMErrorCode SetUpSchurImpl::setUp ( TS  solver)
virtual

Implements SetUpSchur.

Definition at line 1309 of file plastic.cpp.

1309 {
1311 auto simple = mField.getInterface<Simple>();
1312 auto pip_mng = mField.getInterface<PipelineManager>();
1313
1314 SNES snes;
1315 CHKERR TSGetSNES(solver, &snes);
1316 KSP ksp;
1317 CHKERR SNESGetKSP(snes, &ksp);
1318 CHKERR KSPSetFromOptions(ksp);
1319
1320 PC pc;
1321 CHKERR KSPSetFromOptions(ksp);
1322 CHKERR KSPGetPC(ksp, &pc);
1323 PetscBool is_pcfs = PETSC_FALSE;
1324 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1325 if (is_pcfs) {
1326 if (S) {
1329 "Is expected that schur matrix is not allocated. This is "
1330 "possible only is PC is set up twice");
1331 }
1332
1333#ifdef ADD_CONTACT
1334 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1335 A = createDMMatrix(simple->getDM());
1336 P = matDuplicate(A, MAT_DO_NOT_COPY_VALUES);
1337 CHKERR TSSetIJacobian(solver, A, P, TsSetIJacobian, ts_ctx_ptr.get());
1338#endif // ADD_CONTACT
1340 CHKERR MatSetBlockSize(S, SPACE_DIM);
1341
1342 auto set_ops = [&]() {
1344 auto pip_mng = mField.getInterface<PipelineManager>();
1345
1346#ifndef ADD_CONTACT
1347 // Boundary
1348 pip_mng->getOpBoundaryLhsPipeline().push_front(
1349 new OpSchurAssembleBegin());
1350 pip_mng->getOpBoundaryLhsPipeline().push_back(
1352
1353 {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
1354 {SmartPetscObj<Mat>(), S}, {false, false}
1355
1356 ));
1357 // Domain
1358 pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1359 pip_mng->getOpDomainLhsPipeline().push_back(
1361
1362 {"EP", "TAU"}, {nullptr, nullptr}, {SmartPetscObj<AO>(), aoUp},
1363 {SmartPetscObj<Mat>(), S}, {false, false}
1364
1365 ));
1366#else
1367
1368 double eps_stab = 1e-4;
1369 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps_stab", &eps_stab,
1370 PETSC_NULL);
1371
1374 using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
1375
1376 // Boundary
1377 pip_mng->getOpBoundaryLhsPipeline().push_front(
1378 new OpSchurAssembleBegin());
1379 pip_mng->getOpBoundaryLhsPipeline().push_back(
1380 new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
1381 return eps_stab;
1382 }));
1383 pip_mng->getOpBoundaryLhsPipeline().push_back(
1385
1386 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1389 {false, false, false}
1390
1391 ));
1392 // Domain
1393 pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1394 pip_mng->getOpDomainLhsPipeline().push_back(
1396
1397 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr},
1400 {false, false, false}
1401
1402 ));
1403#endif // ADD_CONTACT
1405 };
1406
1407 auto set_assemble_elems = [&]() {
1409 auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
1410 schur_asmb_pre_proc->preProcessHook = [this]() {
1412#ifdef ADD_CONTACT
1413 CHKERR MatZeroEntries(A);
1414 CHKERR MatZeroEntries(P);
1415#endif // ADD_CONTACT
1416 CHKERR MatZeroEntries(S);
1417 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble Begin";
1419 };
1420 auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
1421
1422 schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
1424 MOFEM_LOG("TIMER", Sev::verbose) << "Lhs Assemble End";
1425
1426#ifndef ADD_CONTACT
1428 mField, schur_asmb_post_proc, 1)();
1429#else // ADD_CONTACT
1430 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1431 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1432 // Apply essential constrains to A matrix
1434 mField, schur_asmb_post_proc, 1, A)();
1435 CHKERR MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
1436 CHKERR MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
1437 CHKERR MatAXPY(P, 1, A, SAME_NONZERO_PATTERN);
1438#endif // ADD_CONTACT
1439
1440 // Apply essential constrains to Schur complement
1441 CHKERR MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
1442 CHKERR MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
1444 mField, schur_asmb_post_proc, 1, S, aoUp)();
1445
1447 };
1448 auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
1449 ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
1450 ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
1452 };
1453
1454 auto set_pc = [&]() {
1456 CHKERR PCFieldSplitSetIS(pc, NULL, fieldSplitIS);
1457 CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, S);
1459 };
1460
1461 CHKERR set_ops();
1462 CHKERR set_pc();
1463 CHKERR set_assemble_elems();
1464
1465 } else {
1466 pip_mng->getOpBoundaryLhsPipeline().push_front(new OpSchurAssembleBegin());
1467 pip_mng->getOpBoundaryLhsPipeline().push_back(
1468 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1469 pip_mng->getOpDomainLhsPipeline().push_front(new OpSchurAssembleBegin());
1470 pip_mng->getOpDomainLhsPipeline().push_back(
1471 new OpSchurAssembleEnd<SCHUR_DGESV>({}, {}, {}, {}, {}));
1472 }
1473
1474 // we do not those anymore
1475 subDM.reset();
1476 fieldSplitIS.reset();
1477 // aoUp.reset();
1479}

◆ setUpSubDM()

MoFEMErrorCode SetUpSchurImpl::setUpSubDM ( )
private

Definition at line 953 of file elastic.cpp.

953 {
956 subDM = createDM(mField.get_comm(), "DMMOFEM");
957 CHKERR DMMoFEMCreateSubDM(subDM, simple->getDM(), "SUB");
959 CHKERR DMMoFEMAddElement(subDM, simple->getDomainFEName());
960 auto sub_ents_ptr = boost::make_shared<Range>(subEnts);
961 CHKERR DMMoFEMAddSubFieldRow(subDM, "U", sub_ents_ptr);
962 CHKERR DMMoFEMAddSubFieldCol(subDM, "U", sub_ents_ptr);
963 CHKERR DMSetUp(subDM);
965}

Member Data Documentation

◆ A

SmartPetscObj<Mat> SetUpSchurImpl::A
private
Examples
plastic.cpp.

Definition at line 983 of file contact.cpp.

◆ aoUp

SmartPetscObj< AO > SetUpSchurImpl::aoUp
private
Examples
plastic.cpp.

Definition at line 31 of file SetUpSchurImpl.cpp.

◆ epCorePtr

EshelbianCore* SetUpSchurImpl::epCorePtr
private

Definition at line 36 of file SetUpSchurImpl.cpp.

◆ fieldSplitIS

SmartPetscObj<IS> SetUpSchurImpl::fieldSplitIS
private

IS for split Schur block.

Examples
plastic.cpp.

Definition at line 1305 of file plastic.cpp.

◆ isA00

SmartPetscObj<IS> SetUpSchurImpl::isA00
private

Definition at line 32 of file SetUpSchurImpl.cpp.

◆ isShift

SmartPetscObj<IS> SetUpSchurImpl::isShift
private

Definition at line 33 of file SetUpSchurImpl.cpp.

◆ M

SmartPetscObj<Mat> SetUpSchurImpl::M
private

Definition at line 28 of file SetUpSchurImpl.cpp.

◆ mField

MoFEM::Interface & SetUpSchurImpl::mField
private
Examples
plastic.cpp.

Definition at line 35 of file SetUpSchurImpl.cpp.

◆ P

SmartPetscObj< Mat > SetUpSchurImpl::P
private
Examples
plastic.cpp.

Definition at line 29 of file SetUpSchurImpl.cpp.

◆ S

SmartPetscObj< Mat > SetUpSchurImpl::S
private
Examples
plastic.cpp.

Definition at line 30 of file SetUpSchurImpl.cpp.

◆ subDM

SmartPetscObj< DM > SetUpSchurImpl::subDM
private

field split sub dm

Examples
plastic.cpp.

Definition at line 1304 of file plastic.cpp.

◆ subEnts

Range SetUpSchurImpl::subEnts
private

Definition at line 907 of file elastic.cpp.

◆ volEnts

Range SetUpSchurImpl::volEnts
private

Definition at line 906 of file elastic.cpp.


The documentation for this struct was generated from the following files: