![Logo](MoFEMLogo.png) |
| v0.14.0
|
Go to the documentation of this file.
11 using namespace MoFEM;
13 static char help[] =
"...\n\n";
47 return MatSetValues<AssemblyTypeSelector<PETSC>>(
57 return MatSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
67 return MatSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
79 return MatSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
85 int main(
int argc,
char *argv[]) {
101 DMType dm_name =
"DMMOFEM";
105 auto core_log = logging::core::get();
138 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"VOL",
169 {{simple->getDomainFEName(),
190 auto [mat, block_data_ptr] = shell_data;
193 std::vector<std::string> fields{
"T",
"S",
"O"};
199 {schur_dm, block_dm}, block_data_ptr,
201 fields, {
nullptr,
nullptr,
nullptr}, true
211 using OpMassBlockPreconditionerAssemble =
224 [](
int,
int,
int o) {
return 2 * o; });
225 auto &pip_lhs = pip_mng->getOpDomainLhsPipeline();
227 pip_lhs.push_back(
new OpMassPETSCAssemble(
"V",
"V"));
229 pip_lhs.push_back(
new OpMassPETSCAssemble(
"V",
"T"));
230 pip_lhs.push_back(
new OpMassPETSCAssemble(
"T",
"V"));
231 pip_lhs.push_back(
new OpMassPETSCAssemble(
"S",
"S", close_zero));
232 pip_lhs.push_back(
new OpMassPETSCAssemble(
"S",
"T", beta));
233 pip_lhs.push_back(
new OpMassPETSCAssemble(
"T",
"S", beta));
234 pip_lhs.push_back(
new OpMassPETSCAssemble(
"O",
"O", close_zero));
235 pip_lhs.push_back(
new OpMassPETSCAssemble(
"T",
"O", beta));
236 pip_lhs.push_back(
new OpMassPETSCAssemble(
"O",
"T", beta));
237 pip_lhs.push_back(
new OpMassPETSCAssemble(
"S",
"O", gamma));
238 pip_lhs.push_back(
new OpMassPETSCAssemble(
"O",
"S", gamma));
241 pip_lhs.push_back(
new OpMassBlockAssemble(
"V",
"V"));
243 pip_lhs.push_back(
new OpMassBlockAssemble(
"V",
"T"));
244 pip_lhs.push_back(
new OpMassBlockAssemble(
"T",
"V"));
245 pip_lhs.push_back(
new OpMassBlockAssemble(
"S",
"S", close_zero));
246 pip_lhs.push_back(
new OpMassBlockAssemble(
"S",
"T", beta));
247 pip_lhs.push_back(
new OpMassBlockAssemble(
"T",
"S", beta));
248 pip_lhs.push_back(
new OpMassBlockAssemble(
"O",
"O", close_zero));
249 pip_lhs.push_back(
new OpMassBlockAssemble(
"T",
"O", beta));
250 pip_lhs.push_back(
new OpMassBlockAssemble(
"O",
"T", beta));
251 pip_lhs.push_back(
new OpMassBlockAssemble(
"S",
"O", gamma));
252 pip_lhs.push_back(
new OpMassBlockAssemble(
"O",
"S", gamma));
253 pip_lhs.push_back(
new OpMassBlockPreconditionerAssemble(
"T",
"T"));
255 auto schur_is =
getDMSubData(schur_dm)->getSmartRowIs();
262 {
nullptr,
nullptr,
nullptr},
264 {
nullptr,
nullptr, ao_up}, {
nullptr,
nullptr, S},
266 {
true,
true,
true},
true, block_data_ptr)
273 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
274 MOFEM_LOG(
"Timeline", Sev::inform) <<
"Assemble start";
276 simple->getDomainFEName(),
277 pip_mng->getDomainLhsFE());
278 MOFEM_LOG(
"Timeline", Sev::inform) <<
"Assemble end";
283 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
284 MOFEM_LOG(
"Timeline", Sev::inform) <<
"Mat assemble start";
287 MOFEM_LOG(
"Timeline", Sev::inform) <<
"Mat assemble end";
290 auto get_random_vector = [&](
auto dm) {
293 PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
295 PetscRandomDestroy(&rctx);
298 auto v = get_random_vector(
simple->getDM());
303 auto test = [](
auto msg,
auto y,
double norm0) {
310 CHKERR VecNorm(y, NORM_2, &norm);
314 << msg <<
": norm of difference: " << norm;
315 if (norm >
eps || std::isnan(norm) || std::isinf(norm)) {
316 SETERRQ(PETSC_COMM_WORLD, 1,
"norm of difference is too big");
321 std::vector<int> zero_rows_and_cols = {
325 &*zero_rows_and_cols.begin(), 1, PETSC_NULL,
328 &*zero_rows_and_cols.begin(), 1, PETSC_NULL,
334 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
336 <<
"MatMult(petsc_mat, v, y_petsc) star";
339 <<
"MatMult(petsc_mat, v, y_petsc) end";
342 CHKERR VecNorm(y_petsc, NORM_2, &nrm0);
347 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
349 <<
"MatMult(block_mat, v, y_block) star";
352 <<
"MatMult(block_mat, v, y_block) end";
355 CHKERR VecAXPY(y_petsc, -1.0, y_block);
356 CHKERR test(
"mult", y_petsc, nrm0);
365 CHKERR VecAXPY(y_petsc, -1.0, y_block);
367 CHKERR test(
"mult add", y_petsc, nrm0);
372 CHKERR MatMult(nested_mat,
v, y_nested);
373 CHKERR VecAXPY(y_petsc, -1.0, y_nested);
375 CHKERR test(
"mult nested", y_petsc, nrm0);
378 auto diag_mat = std::get<0>(*nested_data_ptr)[3];
379 auto diag_block_x = get_random_vector(block_dm);
382 CHKERR MatMult(diag_mat, diag_block_x, diag_block_f);
388 CHKERR DMSetMatType(block_dm, MATSHELL);
391 CHKERR DMKSPSetComputeOperators(
393 [](KSP, Mat, Mat,
void *) {
394 MOFEM_LOG(
"WORLD", Sev::inform) <<
"empty operator";
399 auto ksp =
createKSP(m_field.get_comm());
400 CHKERR KSPSetDM(ksp, block_dm);
401 CHKERR KSPSetFromOptions(ksp);
404 auto get_pc = [](
auto ksp) {
406 CHKERR KSPGetPC(ksp, &pc_raw);
412 CHKERR VecZeroEntries(block_solved_x);
413 CHKERR KSPSolve(ksp, diag_block_f, block_solved_x);
416 CHKERR MatMult(diag_mat, block_solved_x, diag_block_f_test);
417 CHKERR VecAXPY(diag_block_f_test, -1.0, diag_block_f);
418 CHKERR test(
"diag solve", diag_block_f_test, nrm0);
420 if (m_field.get_comm_rank() == 0) {
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Data on single entity (This is passed as argument to DataOperator::doWork)
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MPI_Comm & get_comm() const =0
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
@ L2
field with C-1 continuity
UBlasMatrix< double > MatrixDouble
PipelineManager interface.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
auto createKSP(MPI_Comm comm)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Simple interface for fast problem set-up.
auto getDMSubData(DM dm)
Get sub problem data structure.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Deprecated interface functions.
DeprecatedCoreInterface Interface
#define CHKERR
Inline error check.
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, std::vector< SmartPetscObj< AO >> sequence_of_aos, std::vector< SmartPetscObj< Mat >> sequence_of_mats, std::vector< bool > sym_schur, bool symm_op, boost::shared_ptr< BlockStructure > diag_blocks)
Construct a new Op Schur Assemble End object.
auto createDMVector(DM dm)
Get smart vector from DM.
SchurShellMatData createBlockMat(DM dm, boost::shared_ptr< BlockStructure > data)
Create a Schur Mat object.
implementation of Data Operators for Forces and Sources
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Simple interface for fast problem set-up.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
MoFEMErrorCode schurSwitchPreconditioner(boost::shared_ptr< BlockStructure > block_mat_data)
Switch preconditioner.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
constexpr IntegrationType I
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
OpSchurAssembleBase * createOpSchurAssembleBegin()
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< NestSchurData > > createSchurNestedMatrix(boost::shared_ptr< NestSchurData > schur_net_data_ptr)
Create a Mat Diag Blocks object.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Volume finite element base.
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
MoFEMErrorCode DMMoFEMSetBlocMatData(DM dm, boost::shared_ptr< BlockStructure >)
Set data for block mat.
constexpr int SPACE_DIM
[Define dimension]
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
const double v
phase velocity of light in medium (cm/ns)
SmartPetscObj< Mat > block_mat
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define MOFEM_LOG(channel, severity)
Log.
#define CATCH_ERRORS
Catch errors.
EntitiesFieldData::EntData EntData
ForcesAndSourcesCore::UserDataOperator UserDataOperator
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
int main(int argc, char *argv[])
FTensor::Index< 'm', 3 > m
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
MoFEMErrorCode schurSaveBlockMesh(boost::shared_ptr< BlockStructure > block_mat_data, std::string filename)
Save block matrix as a mesh.
SmartPetscObj< Mat > petsc_mat
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
boost::shared_ptr< NestSchurData > getNestSchurData(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.