160 {
161
163#ifndef NDEBUG
164 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble begin -> end";
165#endif
166
168 auto &storage) {
171 for (auto &s : storage) {
172 auto &
m = s.getMat();
174 auto &row_ind = s.getRowInd();
175 auto &col_ind = s.getColInd();
176
177 if (mat_set_values(
M, row_ind.size(), &*row_ind.begin(),
178 col_ind.size(), &*col_ind.begin(),
179 &*
m.data().begin(), ADD_VALUES)) {
181 auto row_ent_it = field_ents->find(s.uidRow);
182 auto col_ent_it = field_ents->find(s.uidCol);
184 if (row_ent_it != field_ents->end())
186 << "Assemble row entity: " << (*row_ent_it)->getName() << " "
187 << (*col_ent_it)->getEntTypeName() << " side "
188 << (*row_ent_it)->getSideNumber();
189 if (col_ent_it != field_ents->end())
191 << "Assemble col entity: " << (*col_ent_it)->getName() << " "
192 << (*col_ent_it)->getEntTypeName() << " side "
193 << (*col_ent_it)->getSideNumber();
195 }
196 }
197 }
198 }
200 };
201
202 auto apply_schur = [&](auto &storage,
203
204 const auto ss,
205
206 auto lo_uid,
auto hi_uid,
double eps) {
208
209 auto add_off_mat = [&](auto row_uid, auto col_uid, auto &row_ind,
212 auto it = storage.template get<SchurL2Mats::uid_mi_tag>().find(
213 boost::make_tuple(row_uid, col_uid));
214 if (it == storage.template get<SchurL2Mats::uid_mi_tag>().end()) {
219 if (idx >= size) {
223 } else {
227 }
229 auto &mat =
p.first->getMat();
230 auto &set_row_ind =
p.first->getRowInd();
231 auto &set_col_ind =
p.first->getColInd();
232#ifndef NDEBUG
233 if (mat.size1() != set_row_ind.size()) {
235 "Wrong size %d != %d", mat.size1(), set_row_ind.size());
236 }
237 if (mat.size2() != set_col_ind.size()) {
239 "Wrong size %d != %d", mat.size2(), set_col_ind.size());
240 }
241#endif
243 noalias(set_row_ind) = row_ind;
244 noalias(set_col_ind) = col_ind;
245 } else {
246 auto &mat = it->getMat();
247#ifndef NDEBUG
250 "Wrong size %d != %d", mat.size1(),
252 }
255 "Wrong size %d != %d", mat.size2(),
257 }
258#endif
260 }
262 };
263
264 auto get_row_view = [&]() {
265 auto row_it =
266 storage.template get<SchurL2Mats::row_mi_tag>().lower_bound(lo_uid);
267 auto hi_row_it =
268 storage.template get<SchurL2Mats::row_mi_tag>().upper_bound(hi_uid);
269 std::vector<const SchurL2Mats *> schur_row_ptr_view;
270 schur_row_ptr_view.reserve(std::distance(row_it, hi_row_it));
271 for (; row_it != hi_row_it; ++row_it) {
272 schur_row_ptr_view.push_back(&*row_it);
273 }
274 return schur_row_ptr_view;
275 };
276
277 auto get_col_view = [&]() {
278 auto col_it =
279 storage.template get<SchurL2Mats::col_mi_tag>().lower_bound(lo_uid);
280 auto hi_col_it =
281 storage.template get<SchurL2Mats::col_mi_tag>().upper_bound(hi_uid);
282 std::vector<const SchurL2Mats *> schur_col_ptr_view;
283 schur_col_ptr_view.reserve(std::distance(col_it, hi_col_it));
284 for (; col_it != hi_col_it; ++col_it) {
285 schur_col_ptr_view.push_back(&*col_it);
286 }
287 return schur_col_ptr_view;
288 };
289
290 auto schur_row_ptr_view = get_row_view();
291 auto schur_col_ptr_view = get_col_view();
292
293 for (auto row_it : schur_row_ptr_view) {
294 if (row_it->uidRow == row_it->uidCol) {
295
297
298 for (auto c_lo : schur_col_ptr_view) {
299
300 auto &row_uid = c_lo->uidRow;
301 if (row_uid == row_it->uidRow)
302 continue;
303
304 auto &cc_off_mat = c_lo->getMat();
306#ifndef NDEBUG
307 if (
invMat.size1() != cc_off_mat.size2()) {
309 "Wrong size %d != %d",
invMat.size1(), cc_off_mat.size2());
310 }
311#endif
313
314 for (auto r_lo : schur_row_ptr_view) {
315 auto &col_uid = r_lo->uidCol;
316
317
318 if (col_uid == row_it->uidRow)
319 continue;
320
321
322 if (
symSchur[ss] && row_uid > col_uid)
323 continue;
324
325 auto &rr_off_mat = r_lo->getMat();
327 rr_off_mat.size2(), false);
328#ifndef NDEBUG
331 "Wrong size %d != %d", rr_off_mat.size1(),
333 }
334#endif
336
337 CHKERR add_off_mat(row_uid, col_uid, c_lo->getRowInd(),
339
340 if (
symSchur[ss] && row_uid != col_uid) {
343 false);
345 CHKERR add_off_mat(col_uid, row_uid, r_lo->getColInd(),
347 }
348 }
349 }
350 }
351 }
352
354 };
355
356 auto erase_factored = [&](auto &storage,
357
358 const auto ss, auto lo_uid, auto hi_uid) {
360
361 auto r_lo =
362 storage.template get<SchurL2Mats::row_mi_tag>().lower_bound(lo_uid);
363 auto r_hi =
364 storage.template get<SchurL2Mats::row_mi_tag>().upper_bound(hi_uid);
365 storage.template get<SchurL2Mats::row_mi_tag>().erase(r_lo, r_hi);
366
367 auto c_lo =
368 storage.template get<SchurL2Mats::col_mi_tag>().lower_bound(lo_uid);
369 auto c_hi =
370 storage.template get<SchurL2Mats::col_mi_tag>().upper_bound(hi_uid);
371 storage.template get<SchurL2Mats::col_mi_tag>().erase(c_lo, c_hi);
372
374 };
375
376 auto assemble_off_diagonal = [&](auto &storage, const auto ss) {
378
379
380 for (
auto &
m : storage) {
382 auto &ind_row =
m.getRowInd();
383 CHKERR AOApplicationToPetsc(ao, ind_row.size(), &*ind_row.begin());
384 auto &ind_col =
m.getColInd();
385 CHKERR AOApplicationToPetsc(ao, ind_col.size(), &*ind_col.begin());
386 }
387 }
388
389
392 }
393
395 };
396
398
399
400
401
402 for (
auto ss = 0; ss !=
fieldsName.size(); ++ss) {
403
404 auto assemble = [&](auto &storage, auto ss, auto lo_uid, auto hi_uid) {
406
407 CHKERR apply_schur(storage, ss, lo_uid, hi_uid,
diagEps[ss]);
408 CHKERR erase_factored(storage, ss, lo_uid, hi_uid);
409 CHKERR assemble_off_diagonal(storage, ss);
411 };
412
415
416 if (row_ents) {
417 for (
auto p = row_ents->pair_begin();
p != row_ents->pair_end(); ++
p) {
418 auto lo_uid =
420 auto hi_uid =
422 CHKERR assemble(storage, ss, lo_uid, hi_uid);
423 }
424 } else {
425
427 field_bit, get_id_for_min_type<MBVERTEX>());
429 field_bit, get_id_for_max_type<MBENTITYSET>());
430 CHKERR assemble(storage, ss, lo_uid, hi_uid);
431 }
432 }
433
434#ifndef NDEBUG
435 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble done";
436#endif
437
439}
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
UBlasMatrix< double > MatrixDouble
UBlasVector< int > VectorInt
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
ForcesAndSourcesCore * getPtrFE() const
MatrixDouble offMatInvDiagOffMat
boost::function< MoFEMErrorCode(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar v[], InsertMode addv)> MatSetValuesRaw
MatSetValuesRaw matSetValuesSchurRaw
MatrixDouble invDiagOffMat
MatrixDouble transOffMatInvDiagOffMat
static boost::ptr_vector< VectorInt > rowIndices
static SchurL2Storage schurL2Storage
static boost::ptr_vector< VectorInt > colIndices
static boost::ptr_vector< MatrixDouble > locMats
Public Types inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator