v0.14.0
Loading...
Searching...
No Matches
Schur.cpp
Go to the documentation of this file.
1/**
2 * @file Schur.cpp
3 * @brief Implementation of Schur Complement
4 * @date 2023-02-03
5 *
6 * @copyright Copyright (c) 2023
7 *
8 */
9
10namespace MoFEM {
11
13boost::ptr_vector<MatrixDouble> SchurL2Mats::locMats;
14boost::ptr_vector<VectorInt> SchurL2Mats::rowIndices;
15boost::ptr_vector<VectorInt> SchurL2Mats::colIndices;
16
17SchurL2Mats::SchurL2Mats(const size_t idx, const UId uid_row, const UId uid_col)
18 : iDX(idx), uidRow(uid_row), uidCol(uid_col) {}
19
22
25 const EntitiesFieldData::EntData &col_data,
26 const MatrixDouble &mat, InsertMode iora) {
27 return MatSetValues<AssemblyTypeSelector<PETSC>>(M, row_data, col_data, mat,
28 iora);
29}
30
33
36 const EntitiesFieldData::EntData &col_data,
37 const MatrixDouble &mat, InsertMode iora) {
39
40 CHKERR matSetValuesPtr(M, row_data, col_data, mat, iora);
41
42 auto get_row_indices = [&]() -> const VectorInt & {
43 if (auto e_ptr = row_data.getFieldEntities()[0]) {
44 if (auto stored_data_ptr =
45 e_ptr->getSharedStoragePtr<EssentialBcStorage>()) {
46 return stored_data_ptr->entityIndices;
47 }
48 }
49 return row_data.getIndices();
50 };
51
52 const auto &row_ind = get_row_indices();
53 const auto &col_ind = col_data.getIndices();
54
55 const auto nb_rows = row_ind.size();
56 const auto nb_cols = col_ind.size();
57
58#ifndef NDEBUG
59 if (mat.size1() != nb_rows) {
60 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
61 "Wrong mat size %d != %d", mat.size1(), nb_rows);
62 }
63 if (mat.size2() != nb_cols) {
64 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
65 "Wrong mat size %d != %d", mat.size2(), nb_cols);
66 }
67#endif // NDEBUG
68
69 const auto idx = SchurL2Mats::schurL2Storage.size();
70 const auto size = SchurL2Mats::locMats.size();
71
72 if (idx >= size) {
73 SchurL2Mats::locMats.push_back(new MatrixDouble());
74 SchurL2Mats::rowIndices.push_back(new VectorInt());
75 SchurL2Mats::colIndices.push_back(new VectorInt());
76 }
77
78 // insert index
79 auto get_uid = [](auto &data) {
80 if (data.getFieldEntities().size() == 1) {
81
82 return data.getFieldEntities()[0]->getLocalUniqueId();
83
84 } else {
85
86 // Is assumed that sum of entities ids gives unique id, that is not true,
87 // but corner case is improbable.
88
89 // @todo: debug version should have test
90
91 auto &uid0 = data.getFieldEntities()[0]->getLocalUniqueId();
92 auto field_id0 = FieldEntity::getFieldBitNumberFromUniqueId(uid0);
93 auto ent0 = FieldEntity::getHandleFromUniqueId(uid0);
94 auto type0 = type_from_handle(ent0);
95 auto id = id_from_handle(ent0);
96
97 for (auto i = 1; i < data.getFieldEntities().size(); ++i) {
98
99 // get entity id from ent
100 id += id_from_handle(
101
102 // get entity handle from unique uid
104 data.getFieldEntities()[i]->getLocalUniqueId())
105
106 );
107 }
108
110 field_id0,
111
112 ent_form_type_and_id(type0, id)
113
114 );
115 }
116 };
117
118 auto uid_row = get_uid(row_data);
119 auto uid_col = get_uid(col_data);
120 auto p = SchurL2Mats::schurL2Storage.emplace(idx, uid_row, uid_col);
121
122 auto get_storage = [&p]() { return const_cast<SchurL2Mats &>(*p.first); };
123
124 if (p.second) {
125
126 auto asmb = [&](auto &sm) {
127 sm.resize(nb_rows, nb_cols, false);
128 noalias(sm) = mat;
129 };
130
131 asmb(get_storage().getMat());
132
133 auto add_indices = [](auto &storage, auto &ind) {
134 storage.resize(ind.size(), false);
135 noalias(storage) = ind;
136 };
137
138 add_indices(get_storage().getRowInd(), row_ind);
139 add_indices(get_storage().getColInd(), col_ind);
140
141 } else {
142
143 auto asmb = [&](auto &sm) {
145
146#ifndef NDEBUG
147 if (sm.size1() != nb_rows) {
148 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
149 "Wrong mat or storage size %d != %d", sm.size1(), nb_rows);
150 }
151 if (sm.size2() != nb_cols) {
152 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
153 "Wrong mat or storage size %d != %d", sm.size2(), nb_cols);
154 }
155#endif // NDEBUG
156
157 switch (iora) {
158 case ADD_VALUES:
159 sm += mat;
160 break;
161 case INSERT_VALUES:
162 noalias(sm) = mat;
163 break;
164 default:
165 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
166 "Assembly type not implemented");
167 }
169 };
170
171 CHKERR asmb(get_storage().getMat());
172 }
173
175}
176
180#ifndef NDEBUG
181 MOFEM_LOG("SELF", Sev::noisy) << "Schur assemble begin";
182#endif
184
185 // auto zero_mats = [&](auto &mats) {
186 // for (auto &m : mats) {
187 // m.resize(0, 0, false);
188 // }
189 // };
190
191 // zero_mats(SchurL2Mats::locMats);
192
194}
195
197 std::vector<std::string> fields_name,
198 std::vector<boost::shared_ptr<Range>> field_ents,
199 std::vector<SmartPetscObj<AO>> sequence_of_aos,
200 std::vector<SmartPetscObj<Mat>> sequence_of_mats,
201 std::vector<bool> sym_schur, std::vector<double> diag_eps, bool symm_op)
202 : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE, symm_op),
203 fieldsName(fields_name), fieldEnts(field_ents),
204 sequenceOfAOs(sequence_of_aos), sequenceOfMats(sequence_of_mats),
205 symSchur(sym_schur), diagEps(diag_eps) {}
206
208 std::vector<std::string> fields_name,
209 std::vector<boost::shared_ptr<Range>> field_ents,
210 std::vector<SmartPetscObj<AO>> sequence_of_aos,
211 std::vector<SmartPetscObj<Mat>> sequence_of_mats,
212 std::vector<bool> sym_schur, bool symm_op)
213 : ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPSPACE, symm_op),
214 fieldsName(fields_name), fieldEnts(field_ents),
215 sequenceOfAOs(sequence_of_aos), sequenceOfMats(sequence_of_mats),
216 symSchur(sym_schur), diagEps(fields_name.size(), 0) {}
217
218template <typename I>
222
224#ifndef NDEBUG
225 MOFEM_LOG("SELF", Sev::noisy) << "Schur assemble begin -> end";
226#endif
227
228 auto assemble = [&](SmartPetscObj<Mat> M, MatSetValuesRaw mat_set_values,
229 auto &storage) {
231 if (M) {
232 for (auto &s : storage) {
233 auto &m = s.getMat();
234 if (m.size1()) {
235 auto &row_ind = s.getRowInd();
236 auto &col_ind = s.getColInd();
237
238 if (mat_set_values(M, row_ind.size(), &*row_ind.begin(),
239 col_ind.size(), &*col_ind.begin(),
240 &*m.data().begin(), ADD_VALUES)) {
241 auto field_ents = getPtrFE()->mField.get_field_ents();
242 auto row_ent_it = field_ents->find(s.uidRow);
243 auto col_ent_it = field_ents->find(s.uidCol);
244 MOFEM_LOG_CHANNEL("SELF");
245 if (row_ent_it != field_ents->end())
246 MOFEM_LOG("SELF", Sev::error)
247 << "Assemble row entity: " << (*row_ent_it)->getName() << " "
248 << (*col_ent_it)->getEntTypeName() << " side "
249 << (*row_ent_it)->getSideNumber();
250 if (col_ent_it != field_ents->end())
251 MOFEM_LOG("SELF", Sev::error)
252 << "Assemble col entity: " << (*col_ent_it)->getName() << " "
253 << (*col_ent_it)->getEntTypeName() << " side "
254 << (*col_ent_it)->getSideNumber();
255 CHK_THROW_MESSAGE(ierr, "MatSetValues");
256 }
257 }
258 }
259 }
261 };
262
263 auto apply_schur = [&](auto &storage,
264
265 const auto ss,
266
267 auto lo_uid, auto hi_uid, double eps) {
269
270 auto add_off_mat = [&](auto row_uid, auto col_uid, auto &row_ind,
271 auto &col_ind, auto &offMatInvDiagOffMat) {
273
274 const auto idx = SchurL2Mats::schurL2Storage.size();
275 const auto size = SchurL2Mats::locMats.size();
276
277 if (idx >= size) {
278 SchurL2Mats::locMats.push_back(new MatrixDouble());
279 SchurL2Mats::rowIndices.push_back(new VectorInt());
280 SchurL2Mats::colIndices.push_back(new VectorInt());
281 }
282
283 auto it = storage.template get<SchurL2Mats::uid_mi_tag>().find(
284 boost::make_tuple(row_uid, col_uid));
285
286 if (it == storage.template get<SchurL2Mats::uid_mi_tag>().end()) {
287
288 auto p = SchurL2Mats::schurL2Storage.emplace(idx, row_uid, col_uid);
289 auto &mat = p.first->getMat();
290 auto &set_row_ind = p.first->getRowInd();
291 auto &set_col_ind = p.first->getColInd();
292
293 set_row_ind.resize(row_ind.size(), false);
294 noalias(set_row_ind) = row_ind;
295 set_col_ind.resize(col_ind.size(), false);
296 noalias(set_col_ind) = col_ind;
297
298 mat.swap(offMatInvDiagOffMat);
299
300 } else {
301 auto &mat = it->getMat();
302#ifndef NDEBUG
303 if (mat.size1() != offMatInvDiagOffMat.size1()) {
304 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
305 "Wrong size %d != %d", mat.size1(),
306 offMatInvDiagOffMat.size1());
307 }
308 if (mat.size2() != offMatInvDiagOffMat.size2()) {
309 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
310 "Wrong size %d != %d", mat.size2(),
311 offMatInvDiagOffMat.size2());
312 }
313#endif // NDEBUG
314 mat += offMatInvDiagOffMat;
315 }
317 };
318
319 auto get_row_view = [&]() {
320 auto row_it =
321 storage.template get<SchurL2Mats::row_mi_tag>().lower_bound(lo_uid);
322 auto hi_row_it =
323 storage.template get<SchurL2Mats::row_mi_tag>().upper_bound(hi_uid);
324 std::vector<const SchurL2Mats *> schur_row_ptr_view;
325 schur_row_ptr_view.reserve(std::distance(row_it, hi_row_it));
326 for (; row_it != hi_row_it; ++row_it) {
327 schur_row_ptr_view.push_back(&*row_it);
328 }
329 return schur_row_ptr_view;
330 };
331
332 auto get_col_view = [&]() {
333 auto col_it =
334 storage.template get<SchurL2Mats::col_mi_tag>().lower_bound(lo_uid);
335 auto hi_col_it =
336 storage.template get<SchurL2Mats::col_mi_tag>().upper_bound(hi_uid);
337 std::vector<const SchurL2Mats *> schur_col_ptr_view;
338 schur_col_ptr_view.reserve(std::distance(col_it, hi_col_it));
339 for (; col_it != hi_col_it; ++col_it) {
340 schur_col_ptr_view.push_back(&*col_it);
341 }
342 return schur_col_ptr_view;
343 };
344
345 auto schur_row_ptr_view = get_row_view();
346 auto schur_col_ptr_view = get_col_view();
347
348 for (auto row_it : schur_row_ptr_view) {
349 if (row_it->uidRow == row_it->uidCol) {
350
351 CHKERR I::invertMat(row_it->getMat(), invMat, eps);
352
353 for (auto c_lo : schur_col_ptr_view) {
354
355 auto &row_uid = c_lo->uidRow;
356 if (row_uid == row_it->uidRow)
357 continue;
358
359 auto &cc_off_mat = c_lo->getMat();
360 invDiagOffMat.resize(cc_off_mat.size1(), invMat.size2(), false);
361#ifndef NDEBUG
362 if (invMat.size1() != cc_off_mat.size2()) {
363 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
364 "Wrong size %d != %d", invMat.size1(), cc_off_mat.size2());
365 }
366#endif // NDEBUG
367
368 // noalias(invDiagOffMat) = prod(cc_off_mat, invMat);
369 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
370 cc_off_mat.size1(), invMat.size2(), cc_off_mat.size2(),
371 1., &*cc_off_mat.data().begin(), cc_off_mat.size2(),
372 &*invMat.data().begin(), invMat.size2(), 0.,
373 &*invDiagOffMat.data().begin(), invDiagOffMat.size2());
374
375 for (auto r_lo : schur_row_ptr_view) {
376 auto &col_uid = r_lo->uidCol;
377
378 // Skip diagonal
379 if (col_uid == row_it->uidRow)
380 continue;
381
382 // If symmetry only run upper off diagonal terms
383 if (symSchur[ss] && row_uid > col_uid)
384 continue;
385
386 auto &rr_off_mat = r_lo->getMat();
387 offMatInvDiagOffMat.resize(invDiagOffMat.size1(),
388 rr_off_mat.size2(), false);
389#ifndef NDEBUG
390 if (rr_off_mat.size1() != invDiagOffMat.size2()) {
391 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
392 "Wrong size %d != %d", rr_off_mat.size1(),
393 invDiagOffMat.size2());
394 }
395#endif // NDEBUG
396
397 // noalias(offMatInvDiagOffMat) = prod(invDiagOffMat, rr_off_mat);
398 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
399 invDiagOffMat.size1(), rr_off_mat.size2(),
400 invDiagOffMat.size2(), 1.,
401 &*invDiagOffMat.data().begin(), invDiagOffMat.size2(),
402 &*rr_off_mat.data().begin(), rr_off_mat.size2(), 0.,
403 &*offMatInvDiagOffMat.data().begin(),
404 offMatInvDiagOffMat.size2());
405
406 CHKERR add_off_mat(row_uid, col_uid, c_lo->getRowInd(),
407 r_lo->getColInd(), offMatInvDiagOffMat);
408
409 if (symSchur[ss] && row_uid != col_uid) {
411 offMatInvDiagOffMat.size1(),
412 false);
414 CHKERR add_off_mat(col_uid, row_uid, r_lo->getColInd(),
415 c_lo->getRowInd(), transOffMatInvDiagOffMat);
416 }
417 }
418 }
419 }
420 }
421
423 };
424
425 auto erase_factored = [&](auto &storage,
426
427 const auto ss, auto lo_uid, auto hi_uid) {
429
430 auto r_lo =
431 storage.template get<SchurL2Mats::row_mi_tag>().lower_bound(lo_uid);
432 auto r_hi =
433 storage.template get<SchurL2Mats::row_mi_tag>().upper_bound(hi_uid);
434 storage.template get<SchurL2Mats::row_mi_tag>().erase(r_lo, r_hi);
435
436 auto c_lo =
437 storage.template get<SchurL2Mats::col_mi_tag>().lower_bound(lo_uid);
438 auto c_hi =
439 storage.template get<SchurL2Mats::col_mi_tag>().upper_bound(hi_uid);
440 storage.template get<SchurL2Mats::col_mi_tag>().erase(c_lo, c_hi);
441
443 };
444
445 auto assemble_off_diagonal = [&](auto &storage, const auto ss) {
447
448 // apply AO
449 for (auto &m : storage) {
450 if (auto ao = sequenceOfAOs[ss]) {
451 auto &ind_row = m.getRowInd();
452 CHKERR AOApplicationToPetsc(ao, ind_row.size(), &*ind_row.begin());
453 auto &ind_col = m.getColInd();
454 CHKERR AOApplicationToPetsc(ao, ind_col.size(), &*ind_col.begin());
455 }
456 }
457
458 // assemble Schur
459 if (sequenceOfMats[ss]) {
460 CHKERR assemble(sequenceOfMats[ss], matSetValuesSchurRaw, storage);
461 }
462
464 };
465
466 auto &storage = SchurL2Mats::schurL2Storage;
467 // Assemble to global matrix
468
469 // Assemble Schur complements
470 // Iterate complements
471 for (auto ss = 0; ss != fieldsName.size(); ++ss) {
472
473 auto assemble = [&](auto &storage, auto ss, auto lo_uid, auto hi_uid) {
475
476 CHKERR apply_schur(storage, ss, lo_uid, hi_uid, diagEps[ss]);
477 CHKERR erase_factored(storage, ss, lo_uid, hi_uid);
478 CHKERR assemble_off_diagonal(storage, ss);
480 };
481
482 auto field_bit = getPtrFE()->mField.get_field_bit_number(fieldsName[ss]);
483 auto row_ents = fieldEnts[ss];
484
485 if (row_ents) {
486 for (auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
487 auto lo_uid =
488 FieldEntity::getLoLocalEntityBitNumber(field_bit, p->first);
489 auto hi_uid =
490 FieldEntity::getHiLocalEntityBitNumber(field_bit, p->second);
491 CHKERR assemble(storage, ss, lo_uid, hi_uid);
492 }
493 } else {
494 // Iterate all entities (typically L2)
496 field_bit, get_id_for_min_type<MBVERTEX>());
498 field_bit, get_id_for_max_type<MBENTITYSET>());
499 CHKERR assemble(storage, ss, lo_uid, hi_uid);
500 }
501 }
502
503#ifndef NDEBUG
504 MOFEM_LOG("SELF", Sev::noisy) << "Schur assemble done";
505#endif
506
508}
509
511 static auto invertMat(MatrixDouble &m, MatrixDouble &inv, double eps) {
513 VectorInt ipiv;
514 VectorDouble lapack_work;
515 const int nb = m.size1();
516
517 if (eps) {
518 for (int cc = 0; cc != nb; ++cc) {
519 m(cc, cc) += eps;
520 }
521 }
522
523 inv.resize(nb, nb, false);
524 inv.clear();
525 auto ptr = &*inv.data().begin();
526 for (int c = 0; c != nb; ++c, ptr += nb + 1)
527 *ptr = -1;
528 ipiv.resize(nb, false);
529 lapack_work.resize(nb * nb, false);
530 const auto info =
531 lapack_dsysv('L', nb, nb, &*m.data().begin(), nb, &*ipiv.begin(),
532 &*inv.data().begin(), nb, &*lapack_work.begin(), nb * nb);
533 if (info != 0)
534 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
535 "Can not invert matrix info = %d", info);
537 };
538};
539
541 static auto invertMat(MatrixDouble &m, MatrixDouble &inv, double eps) {
543 VectorInt ipiv;
544 const auto nb = m.size1();
545#ifndef NDEBUG
546 if (nb != m.size2()) {
547 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
548 "It should be square matrix");
549 }
550#endif
551
552 if (eps) {
553 for (int cc = 0; cc != nb; ++cc) {
554 m(cc, cc) += eps;
555 }
556 }
557
558 inv.resize(nb, nb, false);
559 inv.clear();
560 auto ptr = &*inv.data().begin();
561 for (int c = 0; c != nb; ++c, ptr += nb + 1)
562 *ptr = -1;
563 ipiv.resize(nb, false);
564 const auto info = lapack_dgesv(nb, nb, &*m.data().begin(), nb,
565 &*ipiv.begin(), &*inv.data().begin(), nb);
566 if (info != 0)
567 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
568 "Can not invert matrix info = %d", info);
570 };
571};
572
576 return doWorkImpl<SCHUR_DSYSV>(side, type, data);
577}
578
582 return doWorkImpl<SCHUR_DGESV>(side, type, data);
583}
584
585} // namespace MoFEM
static Index< 'M', 3 > M
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static const double eps
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'm', SPACE_DIM > m
virtual const FieldEntity_multiIndex * get_field_ents() const =0
Get the field ents object.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:176
static __CLPK_integer lapack_dsysv(char uplo, __CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:220
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
uint128_t UId
Unique Id.
Definition: Types.hpp:31
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1869
MoFEMErrorCode schur_mat_set_values_wrap(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:24
auto id_from_handle(const EntityHandle h)
Definition: Templates.hpp:1861
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
Definition: Templates.hpp:1877
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorFieldEntities & getFieldEntities() const
get field entities
const VectorInt & getIndices() const
Get global indices of dofs on entity.
[Storage and set boundary conditions]
static auto getHandleFromUniqueId(const UId uid)
Get the Handle From Unique Id.
static auto getFieldBitNumberFromUniqueId(const UId uid)
Get the Field Bit Number From Unique Id.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: Schur.cpp:177
Assemble Schur complement.
Definition: Schur.hpp:104
MatrixDouble offMatInvDiagOffMat
Definition: Schur.hpp:92
std::vector< boost::shared_ptr< Range > > fieldEnts
Definition: Schur.hpp:84
boost::function< MoFEMErrorCode(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar v[], InsertMode addv)> MatSetValuesRaw
Definition: Schur.hpp:39
static MatSetValuesRaw matSetValuesSchurRaw
Definition: Schur.hpp:41
MoFEMErrorCode doWorkImpl(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: Schur.cpp:220
OpSchurAssembleEndImpl(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=true)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:207
std::vector< bool > symSchur
Definition: Schur.hpp:87
std::vector< std::string > fieldsName
Definition: Schur.hpp:83
MatrixDouble invDiagOffMat
Definition: Schur.hpp:91
std::vector< SmartPetscObj< Mat > > sequenceOfMats
Definition: Schur.hpp:86
std::vector< SmartPetscObj< AO > > sequenceOfAOs
Definition: Schur.hpp:85
MatrixDouble transOffMatInvDiagOffMat
Definition: Schur.hpp:93
std::vector< double > diagEps
Definition: Schur.hpp:88
static auto invertMat(MatrixDouble &m, MatrixDouble &inv, double eps)
Definition: Schur.cpp:541
static auto invertMat(MatrixDouble &m, MatrixDouble &inv, double eps)
Definition: Schur.cpp:511
SchurL2Mats(const size_t idx, const UId uid_row, const UId uid_col)
Definition: Schur.cpp:17
static MatSetValuesPtr matSetValuesPtr
backend assembly function
Definition: Schur.hpp:141
static boost::ptr_vector< VectorInt > rowIndices
Definition: Schur.hpp:182
static SchurL2Storage schurL2Storage
Definition: Schur.hpp:184
multi_index_container< SchurL2Mats, indexed_by< ordered_unique< tag< uid_mi_tag >, composite_key< SchurL2Mats, member< SchurL2Mats, const UId, &SchurL2Mats::uidRow >, member< SchurL2Mats, const UId, &SchurL2Mats::uidCol > > >, ordered_non_unique< tag< row_mi_tag >, member< SchurL2Mats, const UId, &SchurL2Mats::uidRow > >, ordered_non_unique< tag< col_mi_tag >, member< SchurL2Mats, const UId, &SchurL2Mats::uidCol > > > > SchurL2Storage
Definition: Schur.hpp:179
static boost::ptr_vector< VectorInt > colIndices
Definition: Schur.hpp:183
boost::function< MoFEMErrorCode(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)> MatSetValuesPtr
Definition: Schur.hpp:139
auto & getMat() const
Definition: Schur.hpp:132
static MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:35
static boost::ptr_vector< MatrixDouble > locMats
Definition: Schur.hpp:181
auto & getRowInd() const
Definition: Schur.hpp:133
auto & getColInd() const
Definition: Schur.hpp:134
intrusive_ptr for managing petsc objects