MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pmortarassembler.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
13
14#ifdef MFEM_USE_MPI
15
16#include "pmortarassembler.hpp"
17#include "transferutils.hpp"
18
19#include "cut.hpp"
20
21#include "moonolith_bounding_volume_with_span.hpp"
22#include "moonolith_n_tree_mutator_factory.hpp"
23#include "moonolith_n_tree_with_span_mutator_factory.hpp"
24#include "moonolith_n_tree_with_tags_mutator_factory.hpp"
25#include "moonolith_profiler.hpp"
26#include "moonolith_redistribute.hpp"
27#include "moonolith_sparse_matrix.hpp"
28#include "moonolith_tree.hpp"
29#include "par_moonolith.hpp"
30
31#include <memory>
32
33using namespace mfem::internal;
34
35namespace mfem
36{
37
38struct ParMortarAssembler::Impl
39{
40public:
41 MPI_Comm comm;
42 std::shared_ptr<ParFiniteElementSpace> source;
43 std::shared_ptr<ParFiniteElementSpace> destination;
44 std::vector<std::shared_ptr<MortarIntegrator>> integrators;
45
46 std::shared_ptr<HypreParMatrix> coupling_matrix;
47 std::shared_ptr<HypreParMatrix> mass_matrix;
48
49 std::shared_ptr<IterativeSolver> solver;
50
51 bool verbose{false};
52
53 bool assemble_mass_and_coupling_together{true};
54
55 void ensure_solver()
56 {
57 if(!solver) {
58 solver = std::make_shared<BiCGSTABSolver>(comm);
59 solver->SetMaxIter(20000);
60 solver->SetRelTol(1e-6);
61 }
62 }
63
64 bool is_vector_fe() const
65 {
66 bool is_vector_fe = false;
67 for (auto i_ptr : integrators)
68 {
69 if (i_ptr->is_vector_fe())
70 {
71 is_vector_fe = true;
72 break;
73 }
74 }
75
76 return is_vector_fe;
77 }
78
79 BilinearFormIntegrator * new_mass_integrator() const
80 {
81 if (is_vector_fe())
82 {
83 return new VectorFEMassIntegrator();
84 }
85 else
86 {
87 return new MassIntegrator();
88 }
89 }
90};
91
93
94void ParMortarAssembler::SetSolver(const std::shared_ptr<IterativeSolver> &solver)
95{
96 impl_->solver = solver;
97}
98
100{
101 impl_->assemble_mass_and_coupling_together = value;
102}
103
104void ParMortarAssembler::SetMaxSolverIterations(const int max_solver_iterations)
105{
106 impl_->ensure_solver();
107
108 impl_->solver->SetMaxIter(max_solver_iterations);
109}
110
112 const std::shared_ptr<MortarIntegrator> &integrator)
113{
114 impl_->integrators.push_back(integrator);
115}
116
117void ParMortarAssembler::SetVerbose(const bool verbose)
118{
119 impl_->verbose = verbose;
120}
121
122template <int Dimension> class ElementAdapter : public moonolith::Serializable
123{
124public:
125 using Bound = moonolith::AABBWithKDOPSpan<Dimension, double>;
126 using Point = moonolith::Vector<double, Dimension>;
127
128 inline int tag() const { return tag_; }
129
130 const Bound &bound() const { return bound_; }
131
132 Bound &bound() { return bound_; }
133
134 void applyRW(moonolith::Stream &stream)
135 {
136 stream &bound_;
137 stream &element_;
138 stream &element_handle_;
139 }
140
141 ElementAdapter(FiniteElementSpace &fe, const long element,
142 const long element_handle, const int tag)
143 : fe_(&fe), element_(element), element_handle_(element_handle), tag_(tag),
144 dof_map_(nullptr)
145 {
146 assert(element < fe.GetNE());
147
149 fe_->GetMesh()->GetPointMatrix(element, pts);
150
151 Point p;
152 for (int j = 0; j < pts.Width(); ++j)
153 {
154 for (int i = 0; i < pts.Height(); ++i)
155 {
156 p[i] = pts.Elem(i, j);
157 }
158
159 bound_.static_bound() += p;
160 bound_.dynamic_bound() += p;
161 }
162 }
163
164 ElementAdapter()
165 : fe_(nullptr), element_(-1), element_handle_(-1), tag_(-1),
166 dof_map_(nullptr) {}
167
168 inline long handle() const { return element_handle_; }
169
170 inline long element() const { return element_; }
171
172 inline const FiniteElement &get() const
173 {
174 assert(fe_);
175 assert(element_ < fe_->GetNE());
176 return *fe_->GetFE(element_);
177 }
178
179 inline const FiniteElementSpace &space() const
180 {
181 assert(fe_);
182 return *fe_;
183 }
184
185 void set_dof_map(std::vector<long> *ptr) { dof_map_ = ptr; }
186
187 void get_elements_vdofs(Array<int> &vdofs) const
188 {
189 fe_->GetElementVDofs(element_, vdofs);
190
191 if (dof_map_)
192 {
193 assert(dof_map_->size() == vdofs.Size());
194
195 for (int i = 0; i < vdofs.Size(); ++i)
196 {
197 vdofs[i] = dof_map_->at(i);
198 }
199
200 }
201 else
202 {
203 assert(false);
204 }
205 }
206
207private:
209 long element_;
210 long element_handle_;
211 int tag_;
212 Bound bound_;
213 std::vector<long> *dof_map_;
214};
215
216template <int _Dimension> class TreeTraits
217{
218public:
219 enum { Dimension = _Dimension };
220
221 using Bound = moonolith::AABBWithKDOPSpan<Dimension, double>;
222 using DataType = mfem::ElementAdapter<Dimension>;
223};
224
225template <int Dimension>
226class MFEMTree : public moonolith::Tree<TreeTraits<Dimension>>
227{
228public:
229 using Traits = mfem::TreeTraits<Dimension>;
230
231 MFEMTree() {};
232
233 static std::shared_ptr<MFEMTree>
234 New(const int maxElementsXNode = moonolith::DEFAULT_REFINE_MAX_ELEMENTS,
235 const int maxDepth = moonolith::DEFAULT_REFINE_DEPTH)
236 {
237 using namespace moonolith;
238
239 std::shared_ptr<MFEMTree> tree = std::make_shared<MFEMTree>();
240 std::shared_ptr<NTreeWithSpanMutatorFactory<MFEMTree>> factory =
241 std::make_shared<NTreeWithSpanMutatorFactory<MFEMTree>>();
242 factory->set_refine_params(maxElementsXNode, maxDepth);
243 tree->set_mutator_factory(factory);
244 return tree;
245 }
246
247 static std::shared_ptr<MFEMTree>
248 New(const std::shared_ptr<moonolith::Predicate> &predicate,
249 const int maxElementsXNode = moonolith::DEFAULT_REFINE_MAX_ELEMENTS,
250 const int maxDepth = moonolith::DEFAULT_REFINE_DEPTH)
251 {
252 using namespace moonolith;
253
254 if (!predicate)
255 {
256 return New(maxElementsXNode, maxDepth);
257 }
258
259 std::shared_ptr<MFEMTree> tree = std::make_shared<MFEMTree>();
260 std::shared_ptr<NTreeWithTagsMutatorFactory<MFEMTree>> factory =
261 std::make_shared<NTreeWithTagsMutatorFactory<MFEMTree>>(predicate);
262 factory->set_refine_params(maxElementsXNode, maxDepth);
263 tree->set_mutator_factory(factory);
264 return tree;
265 }
266};
267
268class ElementDofMap : public moonolith::Serializable
269{
270public:
271 void read(moonolith::InputStream &is) override
272 {
273 int n;
274 is >> n;
275 global.resize(n);
276 is.read(&global[0], n);
277 }
278
279 void write(moonolith::OutputStream &os) const override
280 {
281 int n = global.size();
282 os << n;
283 os.write(&global[0], n);
284 }
285
286 std::vector<long> global;
287};
288
289class Spaces
290{
291public:
292 explicit Spaces(const moonolith::Communicator &comm) : comm(comm)
293 {
294 must_destroy_attached[0] = false;
295 must_destroy_attached[1] = false;
296 }
297
298 Spaces(const std::shared_ptr<ParFiniteElementSpace> &source,
299 const std::shared_ptr<ParFiniteElementSpace> &destination)
300 {
301 spaces_.reserve(2);
302 spaces_.push_back(source);
303 spaces_.push_back(destination);
304
305 must_destroy_attached[0] = false;
306 must_destroy_attached[1] = false;
307
308 copy_global_dofs(*source, dof_maps_[0]);
309 copy_global_dofs(*destination, dof_maps_[1]);
310 }
311
312 ~Spaces()
313 {
314 Mesh *m = nullptr;
315 FiniteElementCollection *fec = nullptr;
316
317 for (int i = 0; i < spaces_.size(); ++i)
318 {
319 if (spaces_[i] && must_destroy_attached[0])
320 {
321 m = spaces_[i]->GetMesh();
322 fec = const_cast<FiniteElementCollection *>(spaces_[i]->FEColl());
323
324 // make it null
325 spaces_[i] = std::shared_ptr<FiniteElementSpace>();
326
327 delete m;
328 delete fec;
329 }
330 }
331 }
332
333 inline long n_elements() const
334 {
335 long ret = 0;
336 for (auto s : spaces_)
337 {
338 if (s)
339 {
340 ret += s->GetNE();
341 }
342 }
343
344 return ret;
345 }
346
347 inline std::vector<std::shared_ptr<FiniteElementSpace>> &spaces()
348 {
349 return spaces_;
350 }
351
352 inline const std::vector<std::shared_ptr<FiniteElementSpace>> &
353 spaces() const
354 {
355 return spaces_;
356 }
357
358 inline std::vector<ElementDofMap> &dof_map(const int i)
359 {
360 assert(i < 2);
361 assert(i >= 0);
362 return dof_maps_[i];
363 }
364
365 inline const std::vector<ElementDofMap> &dof_map(const int i) const
366 {
367 assert(i < 2);
368 assert(i >= 0);
369 return dof_maps_[i];
370 }
371
372 inline void set_must_destroy_attached(const int index, const bool value)
373 {
374 assert(index < 2);
375 assert(index >= 0);
376 must_destroy_attached[index] = value;
377 }
378
379private:
380 std::vector<std::shared_ptr<FiniteElementSpace>> spaces_;
381 moonolith::Communicator comm;
382 std::vector<ElementDofMap> dof_maps_[2];
383 bool must_destroy_attached[2];
384
385 inline static void copy_global_dofs(ParFiniteElementSpace &fe,
386 std::vector<ElementDofMap> &dof_map)
387 {
388 dof_map.resize(fe.GetNE());
389 Array<int> vdofs;
390 for (int i = 0; i < fe.GetNE(); ++i)
391 {
392 fe.GetElementVDofs(i, vdofs);
393 for (int k = 0; k < vdofs.Size(); ++k)
394 {
395 long g_dof = 0;
396 if (vdofs[k] >= 0)
397 {
398 g_dof = fe.GetGlobalTDofNumber(vdofs[k]);
399 }
400 else
401 {
402 g_dof = -1 - fe.GetGlobalTDofNumber(-1 - vdofs[k]);
403 }
404
405 dof_map[i].global.push_back(g_dof);
406 }
407 }
408 }
409};
410
411template <class Iterator>
412static void write_space(const Iterator &begin, const Iterator &end,
414 const std::vector<ElementDofMap> &dof_map,
415 const int role, moonolith::OutputStream &os)
416{
417 const int dim = space.GetMesh()->Dimension();
418 const long n_elements = std::distance(begin, end);
419
420 std::set<long> nodeIds;
421 std::map<long, long> mapping;
422
423 Array<int> verts;
424 for (Iterator it = begin; it != end; ++it)
425 {
426 const int i = *it;
427 space.GetElementVertices(i, verts);
428
429 for (int j = 0; j < verts.Size(); ++j)
430 {
431 nodeIds.insert(verts[j]);
432 }
433 }
434
435 long n_nodes = nodeIds.size();
436
437 // Estimate for allocation
438 os.request_space((n_elements * 8 + n_nodes * dim) *
439 (sizeof(double) + sizeof(long)));
440
441 auto fe_coll = space.FEColl();
442 const char *name = fe_coll->Name();
443 const int name_lenght = strlen(name);
444 // WRITE 1
445 os << dim << role;
446 os << name_lenght;
447 os.write(name, name_lenght);
448
449 long index = 0;
450 for (auto nodeId : nodeIds)
451 {
452 mapping[nodeId] = index++;
453 }
454
455 // WRITE 2
456 os << n_nodes;
457 // WRITE 6
458 os << n_elements;
459
460 Array<int> vdofs;
461 for (auto node_id : nodeIds)
462 {
463 double *v = space.GetMesh()->GetVertex(node_id);
464 for (int i = 0; i < dim; ++i)
465 {
466 // WRITE 3
467 os << v[i];
468 }
469 }
470
471 for (Iterator it = begin; it != end; ++it)
472 {
473 const int k = *it;
474 space.GetElementVertices(k, verts);
475
476 const int attribute = space.GetAttribute(k);
477 const int e_n_nodes = verts.Size();
478 const int type = space.GetElementType(k);
479 const int order = space.GetOrder(k);
480
481 // WRITE 7
482 os << type << attribute << order << e_n_nodes;
483
484 for (int i = 0; i < e_n_nodes; ++i)
485 {
486 auto it = mapping.find(verts[i]);
487 assert(it != mapping.end());
488
489 int index = it->second;
490
491 // WRITE 8
492 os << index;
493 }
494
495 // WRITE 9
496 os << dof_map.at(k);
497 }
498}
499
500template <class Iterator>
501static void write_element_selection(const Iterator &begin, const Iterator &end,
502 const Spaces &spaces,
503 moonolith::OutputStream &os)
504{
505 if (spaces.spaces().empty())
506 {
507 assert(false);
508 return;
509 }
510
511 auto m = spaces.spaces()[0];
512 std::shared_ptr<FiniteElementSpace> s = nullptr;
513
514 if (spaces.spaces().size() > 1)
515 {
516 s = spaces.spaces()[1];
517 }
518
519 std::vector<long> source_selection;
520 std::vector<long> destination_selection;
521
522 bool met_destination_selection = false;
523
524 for (Iterator it = begin; it != end; ++it)
525 {
526 long index = *it;
527
528 if (m && index >= m->GetNE())
529 {
530 index -= m->GetNE();
531 destination_selection.push_back(index);
532 }
533 else if (!m)
534 {
535 met_destination_selection = true;
536 destination_selection.push_back(index);
537 }
538 else
539 {
540 assert(!met_destination_selection);
541 assert(index < m->GetNE());
542 source_selection.push_back(index);
543 }
544 }
545
546 const bool has_source = !source_selection.empty();
547 const bool has_destination = !destination_selection.empty();
548
549 os << has_source << has_destination;
550
551 if (has_source)
552 {
553 write_space(source_selection.begin(), source_selection.end(), *m,
554 spaces.dof_map(0), 0, os);
555 }
556
557 if (has_destination)
558 {
559 write_space(destination_selection.begin(), destination_selection.end(), *s,
560 spaces.dof_map(1), 1, os);
561 }
562}
563
564static FiniteElementCollection *FECollFromName(const std::string &comp_name)
565{
566 return FiniteElementCollection::New(comp_name.c_str());
567}
568
569static void read_space(moonolith::InputStream &is,
570 std::shared_ptr<FiniteElementSpace> &space,
571 std::vector<ElementDofMap> &dof_map)
572{
573 using namespace std;
574
575 // READ 1
576 int dim, role, name_lenght;
577 is >> dim >> role;
578 is >> name_lenght;
579
580 std::string name(name_lenght, 0);
581 is.read(&name[0], name_lenght);
582
583 // READ 2
584 long n_nodes;
585 is >> n_nodes;
586
587 // READ 6
588 long n_elements;
589 is >> n_elements;
590
591 auto fe_coll = FECollFromName(name);
592 auto mesh_ptr = new Mesh(dim, n_nodes, n_elements);
593
594 for (long i = 0; i < n_nodes; ++i)
595 {
596 double v[3];
597 for (int i = 0; i < dim; ++i)
598 {
599 // READ 3
600 is >> v[i];
601 }
602
603 mesh_ptr->AddVertex(v);
604 }
605
606 dof_map.resize(n_elements);
607 std::vector<int> e2v;
608 for (long i = 0; i < n_elements; ++i)
609 {
610 // READ 7
611 int type, attribute, order, e_n_nodes;
612 is >> type >> attribute >> order >> e_n_nodes;
613 e2v.resize(e_n_nodes);
614 int index;
615 for (int i = 0; i < e_n_nodes; ++i)
616 {
617 // READ 8
618 is >> index;
619 e2v[i] = index;
620 }
621
622 mesh_ptr->AddElement(NewElem(type, &e2v[0], attribute));
623 // READ 9
624 is >> dof_map.at(i);
625 }
626
627 // if(mesh_ptr->Dimension() == 3) {
628 Finalize(*mesh_ptr, true);
629 // }
630
631 space = make_shared<FiniteElementSpace>(mesh_ptr, fe_coll);
632}
633
634static void read_spaces(moonolith::InputStream &is, Spaces &spaces)
635{
636 bool has_source, has_destination;
637 is >> has_source >> has_destination;
638
639 spaces.spaces().resize(2);
640
641 if (has_source)
642 {
643 read_space(is, spaces.spaces()[0], spaces.dof_map(0));
644 spaces.set_must_destroy_attached(0, true);
645 }
646 else
647 {
648 spaces.spaces()[0] = nullptr;
649 }
650
651 if (has_destination)
652 {
653 read_space(is, spaces.spaces()[1], spaces.dof_map(1));
654 spaces.set_must_destroy_attached(1, true);
655 }
656 else
657 {
658 spaces.spaces()[1] = nullptr;
659 }
660}
661
662template <int Dimensions, class Fun>
663static bool Assemble(moonolith::Communicator &comm,
664 std::shared_ptr<ParFiniteElementSpace> &source,
665 std::shared_ptr<ParFiniteElementSpace> &destination,
666 Fun process_fun, const moonolith::SearchSettings &settings,
667 const bool verbose)
668{
669 using namespace moonolith;
670
671 typedef mfem::MFEMTree<Dimensions> NTreeT;
672 typedef typename NTreeT::DataContainer DataContainer;
673 typedef typename NTreeT::DataType Adapter;
674
675 long maxNElements = settings.max_elements;
676 long maxDepth = settings.max_depth;
677
678 const int n_elements_source = source->GetNE();
679 const int n_elements_destination = destination->GetNE();
680 const int n_elements = n_elements_source + n_elements_destination;
681
682 auto predicate = std::make_shared<MasterAndSlave>();
683 predicate->add(0, 1);
684
685 MOONOLITH_EVENT_BEGIN("create_adapters");
686
687 std::shared_ptr<NTreeT> tree = NTreeT::New(predicate, maxNElements, maxDepth);
688 tree->reserve(n_elements);
689
690 std::shared_ptr<Spaces> local_spaces =
691 std::make_shared<Spaces>(source, destination);
692
693 int offset = 0;
694 int space_num = 0;
695
696 for (auto s : local_spaces->spaces())
697 {
698 if (s)
699 {
700 for (int i = 0; i < s->GetNE(); ++i)
701 {
702 Adapter a(*s, i, offset + i, space_num);
703 a.set_dof_map(&local_spaces->dof_map(space_num)[i].global);
704 tree->insert(a);
705 }
706
707 offset += s->GetNE();
708 }
709
710 ++space_num;
711 }
712
713 tree->root()->bound().static_bound().enlarge(1e-6);
714
715 MOONOLITH_EVENT_END("create_adapters");
716
717 // Just to have an indexed-storage
718 std::map<long, std::shared_ptr<Spaces>> spaces;
719 std::map<long, std::vector<std::shared_ptr<Spaces>>> migrated_spaces;
720
721 auto read = [&spaces, &migrated_spaces,
722 comm](const long ownerrank, const long /*senderrank*/,
723 bool is_forwarding, DataContainer &data, InputStream &in)
724 {
725 CHECK_STREAM_READ_BEGIN("vol_proj", in);
726
727 std::shared_ptr<Spaces> proc_space = std::make_shared<Spaces>(comm);
728
729 read_spaces(in, *proc_space);
730
731 if (!is_forwarding)
732 {
733 assert(!spaces[ownerrank]);
734 spaces[ownerrank] = proc_space;
735 }
736 else
737 {
738 migrated_spaces[ownerrank].push_back(proc_space);
739 }
740
741 data.reserve(data.size() + proc_space->n_elements());
742
743 int space_num = 0;
744 long offset = 0;
745 for (auto s : proc_space->spaces())
746 {
747 if (s)
748 {
749 for (int i = 0; i < s->GetNE(); ++i)
750 {
751 data.push_back(Adapter(*s, i, offset + i, space_num));
752 data.back().set_dof_map(&proc_space->dof_map(space_num)[i].global);
753 }
754
755 offset += s->GetNE();
756 }
757
758 ++space_num;
759 }
760
761 CHECK_STREAM_READ_END("vol_proj", in);
762 };
763
764 auto write = [&local_spaces, &spaces,
765 &comm](const long ownerrank, const long /*recvrank*/,
766 const std::vector<long>::const_iterator &begin,
767 const std::vector<long>::const_iterator &end,
768 const DataContainer &/*data*/, OutputStream &out)
769 {
770 CHECK_STREAM_WRITE_BEGIN("vol_proj", out);
771
772 if (ownerrank == comm.rank())
773 {
774 write_element_selection(begin, end, *local_spaces, out);
775 }
776 else
777 {
778 auto it = spaces.find(ownerrank);
779 assert(it != spaces.end());
780 std::shared_ptr<Spaces> spaceptr = it->second;
781 assert(std::distance(begin, end) > 0);
782 write_element_selection(begin, end, *spaceptr, out);
783 }
784
785 CHECK_STREAM_WRITE_END("vol_proj", out);
786 };
787
788 long n_false_positives = 0, n_intersections = 0;
789 auto fun = [&n_false_positives, &n_intersections,
790 &process_fun](Adapter &source, Adapter &destination) -> bool
791 {
792 bool ok = process_fun(source, destination);
793
794 if (ok)
795 {
796 n_intersections++;
797 return true;
798 }
799 else
800 {
801 n_false_positives++;
802 return false;
803 }
804
805 return true;
806 };
807
808 moonolith::search_and_compute(comm, tree, predicate, read, write, fun,
809 settings);
810
811 if (verbose)
812 {
813 long n_total_candidates = n_intersections + n_false_positives;
814
815 long n_collection[3] = {n_intersections, n_total_candidates,
816 n_false_positives
817 };
818 comm.all_reduce(n_collection, 3, moonolith::MPISum());
819
820 if (comm.is_root())
821 {
822 mfem::out << "n_intersections: " << n_collection[0]
823 << ", n_total_candidates: " << n_collection[1]
824 << ", n_false_positives: " << n_collection[2] << std::endl;
825 }
826 }
827
828 return true;
829}
830
831static void add_matrix(const Array<int> &destination_vdofs, const Array<int> &source_vdofs,
832 const DenseMatrix &elem_mat, moonolith::SparseMatrix<double> &mat_buffer)
833{
834 for (int i = 0; i < destination_vdofs.Size(); ++i)
835 {
836 long dof_I = destination_vdofs[i];
837
838 double sign_I = 1.0;
839
840 if (dof_I < 0)
841 {
842 sign_I = -1.0;
843 dof_I = -dof_I - 1;
844 }
845
846 for (int j = 0; j < source_vdofs.Size(); ++j)
847 {
848 long dof_J = source_vdofs[j];
849
850 double sign_J = sign_I;
851
852 if (dof_J < 0)
853 {
854 sign_J = -sign_I;
855 dof_J = -dof_J - 1;
856 }
857
858 mat_buffer.add(dof_I, dof_J, sign_J * elem_mat.Elem(i, j));
859 }
860 }
861}
862
863std::shared_ptr<HypreParMatrix> convert_to_hypre_matrix(
864 const std::vector<moonolith::Integer> &destination_ranges,
865 HYPRE_BigInt *s_offsets,
866 HYPRE_BigInt *m_offsets,
867 moonolith::SparseMatrix<double> &mat_buffer) {
868
869 auto && comm = mat_buffer.comm();
870
871 moonolith::Redistribute<moonolith::SparseMatrix<double>> redist(comm.get());
872 redist.apply(destination_ranges, mat_buffer, moonolith::AddAssign<double>());
873
874 std::vector<int> I(s_offsets[1] - s_offsets[0] + 1);
875 I[0] = 0;
876
877 std::vector<HYPRE_Int> J;
878 std::vector<double> data;
879 J.reserve(mat_buffer.n_local_entries());
880 data.reserve(J.size());
881
882 for (auto it = mat_buffer.iter(); it; ++it)
883 {
884 I[it.row() - s_offsets[0] + 1]++;
885 J.push_back(it.col());
886 data.push_back(*it);
887 }
888
889 for (int i = 1; i < I.size(); ++i)
890 {
891 I[i] += I[i - 1];
892 }
893
894 return std::make_shared<HypreParMatrix>(
895 comm.get(), s_offsets[1] - s_offsets[0], mat_buffer.rows(),
896 mat_buffer.cols(), &I[0], &J[0], &data[0], s_offsets, m_offsets);
897}
898
899
900int order_multiplier(const Geometry::Type type, const int dim)
901{
902 return
903 (type == Geometry::TRIANGLE || type == Geometry::TETRAHEDRON ||
904 type == Geometry::SEGMENT)? 1 : dim;
905}
906
907template <int Dimensions>
908static bool
909Assemble(moonolith::Communicator &comm,
910 ParMortarAssembler::Impl &impl,
911 const moonolith::SearchSettings &settings, const bool verbose)
912{
913 auto &source = impl.source;
914 auto &destination = impl.destination;
915
916 auto &integrators = impl.integrators;
917 auto &pmat = impl.coupling_matrix;
918
919
920 bool lump_mass = false;
921 int max_q_order = 0;
922
923 for (auto i_ptr : integrators)
924 {
925 max_q_order = std::max(i_ptr->GetQuadratureOrder(), max_q_order);
926 }
927
928 const int dim = source->GetMesh()->Dimension();
929 std::shared_ptr<Cut> cut = NewCut(dim);
930 if (!cut)
931 {
932 assert(false && "NOT Supported!");
933 return false;
934 }
935
936 Array<int> source_vdofs, destination_vdofs;
937 DenseMatrix elemmat;
938 DenseMatrix cumulative_elemmat;
939 IntegrationRule src_ir;
940 IntegrationRule dest_ir;
941
942 double local_element_matrices_sum = 0.0;
943
944 const auto m_global_n_dofs = source->GlobalVSize();
945 auto *m_offsets = source->GetDofOffsets();
946
947 const auto s_global_n_dofs = destination->GlobalVSize();
948 auto *s_offsets = destination->GetDofOffsets();
949
950 moonolith::SparseMatrix<double> mat_buffer(comm);
951 mat_buffer.set_size(s_global_n_dofs, m_global_n_dofs);
952
953 moonolith::SparseMatrix<double> mass_mat_buffer(comm);
954 mass_mat_buffer.set_size(s_global_n_dofs, s_global_n_dofs);
955
956 std::unique_ptr<BilinearFormIntegrator> mass_integr(impl.new_mass_integrator());
957
958 auto fun = [&](const ElementAdapter<Dimensions> &source,
959 const ElementAdapter<Dimensions> &destination) -> bool
960 {
961 const auto &src = source.space();
962 const auto &dest = destination.space();
963
964 const int src_index = source.element();
965 const int dest_index = destination.element();
966
967 auto &src_fe = *src.GetFE(src_index);
968 auto &dest_fe = *dest.GetFE(dest_index);
969
970 ElementTransformation &dest_Trans =
971 *dest.GetElementTransformation(dest_index);
972
973 // Quadrature order mangling
974
975 int src_order_mult = order_multiplier(src_fe.GetGeomType(), Dimensions);
976 int dest_order_mult = order_multiplier(dest_fe.GetGeomType(), Dimensions);
977
978 const int src_order = src_order_mult * src_fe.GetOrder();
979 const int dest_order = dest_order_mult * dest_fe.GetOrder();
980
981 int contraction_order = src_order + dest_order;
982 if(impl.assemble_mass_and_coupling_together) {
983 contraction_order = std::max(contraction_order, 2 * dest_order);
984 }
985
986 const int order = contraction_order + dest_order_mult * dest_Trans.OrderW() + max_q_order;
987
988 cut->SetIntegrationOrder(order);
989
990 if (cut->BuildQuadrature(src, src_index, dest, dest_index, src_ir,
991 dest_ir))
992 {
993 // make reference quadratures
994 ElementTransformation &src_Trans =
995 *src.GetElementTransformation(src_index);
996
997 source.get_elements_vdofs(source_vdofs);
998 destination.get_elements_vdofs(destination_vdofs);
999
1000 bool first = true;
1001 for (auto i_ptr : integrators)
1002 {
1003 if (first)
1004 {
1005 i_ptr->AssembleElementMatrix(src_fe, src_ir, src_Trans, dest_fe,
1006 dest_ir, dest_Trans, cumulative_elemmat);
1007 first = false;
1008 }
1009 else
1010 {
1011 i_ptr->AssembleElementMatrix(src_fe, src_ir, src_Trans, dest_fe,
1012 dest_ir, dest_Trans, elemmat);
1013 cumulative_elemmat += elemmat;
1014 }
1015 }
1016
1017 local_element_matrices_sum += Sum(cumulative_elemmat);
1018 add_matrix(destination_vdofs, source_vdofs, cumulative_elemmat, mat_buffer);
1019
1020 if(impl.assemble_mass_and_coupling_together) {
1021 mass_integr->SetIntRule(&dest_ir);
1022 mass_integr->AssembleElementMatrix(dest_fe, dest_Trans, elemmat);
1023
1024 if(lump_mass) {
1025 int n = destination_vdofs.Size();
1026
1027 for(int i = 0; i < n; ++i) {
1028 double row_sum = 0.;
1029 for(int j = 0; j < n; ++j) {
1030 row_sum += elemmat(i, j);
1031 elemmat(i, j) = 0.;
1032 }
1033
1034 elemmat(i, i) = row_sum;
1035 }
1036 }
1037
1038 add_matrix(destination_vdofs, destination_vdofs, elemmat, mass_mat_buffer);
1039 }
1040
1041 return true;
1042 }
1043 else
1044 {
1045 return false;
1046 }
1047 };
1048
1049 if (!Assemble<Dimensions>(comm, source, destination, fun, settings,
1050 verbose))
1051 {
1052 return false;
1053 }
1054
1055 if (verbose)
1056 {
1057 double volumes[2] = {local_element_matrices_sum};
1058 comm.all_reduce(volumes, 2, moonolith::MPISum());
1059
1060 cut->Describe();
1061
1062 if (comm.is_root())
1063 {
1064 mfem::out << "sum(B): " << volumes[0] << std::endl;
1065 }
1066 }
1067
1068 std::vector<moonolith::Integer> destination_ranges(comm.size() + 1, 0);
1069
1070 std::copy(s_offsets, s_offsets + 2, destination_ranges.begin() + comm.rank());
1071 comm.all_reduce(&destination_ranges[0], destination_ranges.size(),
1072 moonolith::MPIMax());
1073
1075 destination_ranges,
1076 s_offsets,
1077 m_offsets,
1078 mat_buffer);
1079
1080 if(impl.assemble_mass_and_coupling_together) {
1081 impl.mass_matrix = convert_to_hypre_matrix(destination_ranges, s_offsets, s_offsets, mass_mat_buffer);
1082 }
1083
1084 return true;
1085}
1086
1088 const std::shared_ptr<ParFiniteElementSpace> &source,
1089 const std::shared_ptr<ParFiniteElementSpace> &destination)
1090 : impl_(new Impl())
1091{
1092 impl_->comm = source->GetComm();
1093 impl_->source = source;
1094 impl_->destination = destination;
1095}
1096
1097bool ParMortarAssembler::Assemble(std::shared_ptr<HypreParMatrix> &pmat)
1098{
1099 assert(!impl_->integrators.empty() &&
1100 "it must have at least on integrator see class MortarIntegrator");
1101
1102 moonolith::SearchSettings settings;
1103 // settings.set("disable_redistribution", moonolith::Boolean(true));
1104 // settings.set("disable_asynch", moonolith::Boolean(true));
1105
1106 bool ok = false;
1107
1108 moonolith::Communicator comm(impl_->comm);
1109 if (impl_->source->GetMesh()->Dimension() == 2)
1110 {
1111 ok = mfem::Assemble<2>(comm, *impl_, settings,
1112 impl_->verbose);
1113 }
1114
1115 if (impl_->source->GetMesh()->Dimension() == 3)
1116 {
1117 ok = mfem::Assemble<3>(comm, *impl_, settings,
1118 impl_->verbose);
1119 }
1120
1121 pmat = impl_->coupling_matrix;
1122 return ok;
1123}
1124
1126 ParGridFunction &dest_fun)
1127{
1128
1129 return Update() && Apply(src_fun, dest_fun);
1130}
1131
1133 ParGridFunction &dest_fun)
1134{
1135 const bool verbose = impl_->verbose;
1136
1137 if (!impl_->coupling_matrix)
1138 {
1139 if (!Update())
1140 {
1141 return false;
1142 }
1143 }
1144
1145 auto &B = *impl_->coupling_matrix;
1146 auto &D = *impl_->mass_matrix;
1147
1148 auto &P_source = *impl_->source->Dof_TrueDof_Matrix();
1149 auto &P_destination = *impl_->destination->Dof_TrueDof_Matrix();
1150
1151 impl_->ensure_solver();
1152
1153 // BlockILU prec(D);
1154 // impl_->solver->SetPreconditioner(prec);
1155
1156 impl_->solver->SetOperator(D);
1157
1158 if(verbose) {
1159 impl_->solver->SetPrintLevel(3);
1160 }
1161
1162 Vector P_x_src_fun(B.Width());
1163 P_source.MultTranspose(src_fun, P_x_src_fun);
1164
1165 Vector B_x_src_fun(B.Height());
1166 B_x_src_fun = 0.0;
1167
1168 B.Mult(P_x_src_fun, B_x_src_fun);
1169
1170 Vector R_x_dest_fun(D.Height());
1171 R_x_dest_fun = 0.0;
1172
1173 impl_->solver->Mult(B_x_src_fun, R_x_dest_fun);
1174
1175 P_destination.Mult(R_x_dest_fun, dest_fun);
1176
1177 dest_fun.Update();
1178
1179 // Sanity check!
1180 int converged = impl_->solver->GetFinalNorm() < 1e-5;
1181 MPI_Allreduce(MPI_IN_PLACE, &converged, 1, MPI_INT, MPI_MIN, impl_->comm);
1182 return converged;
1183}
1184
1186{
1187 using namespace std;
1188 const bool verbose = impl_->verbose;
1189 static const bool dof_transformation = true;
1190
1191 moonolith::Communicator comm(impl_->comm);
1192
1193 if (verbose)
1194 {
1195
1196 moonolith::root_describe(
1197 "--------------------------------------------------------"
1198 "\nAssembly begin: ",
1199 comm, mfem::out);
1200 }
1201
1202 if (!Assemble(impl_->coupling_matrix))
1203 {
1204 return false;
1205 }
1206
1207 if (verbose)
1208 {
1209 moonolith::root_describe(
1210 "\nAssembly end: "
1211 "--------------------------------------------------------",
1212 comm, mfem::out);
1213 }
1214
1215 if (dof_transformation)
1216 {
1217 impl_->coupling_matrix.reset(RAP(impl_->destination->Dof_TrueDof_Matrix(),
1218 impl_->coupling_matrix.get(),
1219 impl_->source->Dof_TrueDof_Matrix()));
1220 }
1221
1222 auto &B = *impl_->coupling_matrix;
1223
1224 if(!impl_->assemble_mass_and_coupling_together) {
1225 ParBilinearForm b_form(impl_->destination.get());
1226 b_form.AddDomainIntegrator(impl_->new_mass_integrator());
1227 b_form.Assemble();
1228 b_form.Finalize();
1229 impl_->mass_matrix = shared_ptr<HypreParMatrix>(b_form.ParallelAssemble());
1230 }
1231
1232
1233 comm.barrier();
1234
1235 if (verbose && comm.is_root())
1236 {
1237 mfem::out << "P in R^(" << impl_->source->Dof_TrueDof_Matrix()->Height();
1238 mfem::out << " x " << impl_->source->Dof_TrueDof_Matrix()->Width() << ")\n";
1239 mfem::out << "Q in R^("
1240 << impl_->destination->Dof_TrueDof_Matrix()->Height();
1241 mfem::out << " x " << impl_->destination->Dof_TrueDof_Matrix()->Width()
1242 << ")\n";
1243 }
1244
1245 if (dof_transformation && impl_->assemble_mass_and_coupling_together)
1246 {
1247 impl_->mass_matrix.reset(RAP(impl_->destination->Dof_TrueDof_Matrix(),
1248 impl_->mass_matrix.get(),
1249 impl_->destination->Dof_TrueDof_Matrix()));
1250 }
1251
1252 auto &D = *impl_->mass_matrix;
1253
1254 comm.barrier();
1255
1256 if (verbose && comm.is_root())
1257 {
1258 mfem::out << "--------------------------------------------------------"
1259 << std::endl;
1260 mfem::out << "B in R^(" << B.GetGlobalNumRows() << " x "
1261 << B.GetGlobalNumCols() << ")" << std::endl;
1262 mfem::out << "D in R^(" << D.GetGlobalNumRows() << " x "
1263 << D.GetGlobalNumCols() << ")" << std::endl;
1264 mfem::out << "--------------------------------------------------------"
1265 << std::endl;
1266 }
1267
1268 comm.barrier();
1269
1270 if (verbose)
1271 {
1272 Vector v(D.Width());
1273 v = 1.0;
1274
1275 Vector Dv(D.Height());
1276 D.Mult(v, Dv);
1277
1278 double sum_Dv = Dv.Sum();
1279 comm.all_reduce(&sum_Dv, 1, MPI_SUM);
1280
1281 if (comm.is_root())
1282 {
1283 mfem::out << "sum(D): " << sum_Dv << std::endl;
1284 }
1285 }
1286
1287 return true;
1288}
1289
1290} // namespace mfem
1291
1292#endif // MFEM_USE_MPI
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Abstract base class BilinearFormIntegrator.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
virtual real_t & Elem(int i, int j)
Returns reference to a_{ij}.
Definition densemat.cpp:210
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name.
Definition fe_coll.cpp:126
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:287
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:740
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
Abstract class for all finite elements.
Definition fe_base.hpp:239
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
Mesh data type.
Definition mesh.hpp:56
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
void GetPointMatrix(int i, DenseMatrix &pointmat) const
Definition mesh.cpp:7326
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GetGlobalTDofNumber(int ldof) const
Returns the global tdof number of the given local degree of freedom.
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Update() override
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition pgridfunc.cpp:90
bool Apply(const ParGridFunction &src_fun, ParGridFunction &dest_fun)
transfer a function from source to destination. It requires that the Update function is called before
void SetSolver(const std::shared_ptr< IterativeSolver > &solver)
Changes the solver to be used for solving the mass-matrix linear system.
void SetMaxSolverIterations(const int max_solver_iterations)
Control the maximum numbers of conjugate gradients steps for mass matrix inversion.
bool Assemble(std::shared_ptr< HypreParMatrix > &B)
assembles the coupling matrix B. B : source -> destination If u is a coefficient associated with sour...
ParMortarAssembler(const std::shared_ptr< ParFiniteElementSpace > &source, const std::shared_ptr< ParFiniteElementSpace > &destination)
constructs the object with source and destination spaces
void SetAssembleMassAndCouplingTogether(const bool value)
Control if the Mass matrix is computed together with the coupling operator every time.
bool Transfer(const ParGridFunction &src_fun, ParGridFunction &dest_fun)
transfer a function from source to destination. if the transfer is to be performed multiple times use...
void SetVerbose(const bool verbose)
Expose process details with verbose output.
void AddMortarIntegrator(const std::shared_ptr< MortarIntegrator > &integrator)
This method must be called before Assemble or Transfer. It will assemble the operator in all intersec...
bool Update()
assembles the various components necessary for the transfer. To be called before calling the Apply fu...
Data type point element.
Definition point.hpp:23
Vector data type.
Definition vector.hpp:80
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
int dim
Definition ex24.cpp:53
void source(const Vector &x, Vector &f)
Definition ex25.cpp:620
HYPRE_Int HYPRE_BigInt
int index(int i, int j, int nx, int ny)
Definition life.cpp:235
real_t a
Definition lissajous.cpp:41
string space
void write(std::ostream &os, T value)
Write 'value' to stream.
Definition binaryio.hpp:30
T read(std::istream &is)
Read a value from the stream and return it.
Definition binaryio.hpp:37
std::shared_ptr< Cut > NewCut(const int dim)
Definition cut.cpp:436
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
int order_multiplier(const Geometry::Type type, const int dim)
std::shared_ptr< HypreParMatrix > convert_to_hypre_matrix(const std::vector< moonolith::Integer > &destination_ranges, HYPRE_BigInt *s_offsets, HYPRE_BigInt *m_offsets, moonolith::SparseMatrix< double > &mat_buffer)
real_t p(const Vector &x, real_t t)
RefCoord s[3]
void pts(int iphi, int t, real_t x[])