MFEM  v4.5.2
Finite element discretization library
pmortarassembler.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 
12 #include "../../config/config.hpp"
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 
33 using namespace mfem::internal;
34 
35 namespace mfem
36 {
37 
38 struct ParMortarAssembler::Impl
39 {
40 public:
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 
94 void 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 
104 void 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 
117 void ParMortarAssembler::SetVerbose(const bool verbose)
118 {
119  impl_->verbose = verbose;
120 }
121 
122 template <int Dimension> class ElementAdapter : public moonolith::Serializable
123 {
124 public:
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 
148  DenseMatrix pts;
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 
207 private:
208  FiniteElementSpace *fe_;
209  long element_;
210  long element_handle_;
211  int tag_;
212  Bound bound_;
213  std::vector<long> *dof_map_;
214 };
215 
216 template <int _Dimension> class TreeTraits
217 {
218 public:
219  enum { Dimension = _Dimension };
220 
221  using Bound = moonolith::AABBWithKDOPSpan<Dimension, double>;
222  using DataType = mfem::ElementAdapter<Dimension>;
223 };
224 
225 template <int Dimension>
226 class MFEMTree : public moonolith::Tree<TreeTraits<Dimension>>
227 {
228 public:
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 
268 class ElementDofMap : public moonolith::Serializable
269 {
270 public:
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 
289 class Spaces
290 {
291 public:
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 
379 private:
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 
411 template <class Iterator>
412 static void write_space(const Iterator &begin, const Iterator &end,
413  FiniteElementSpace &space,
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 
500 template <class Iterator>
501 static 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 
564 static FiniteElementCollection *FECollFromName(const std::string &comp_name)
565 {
566  return FiniteElementCollection::New(comp_name.c_str());
567 }
568 
569 static 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 
634 static 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 
662 template <int Dimensions, class Fun>
663 static 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 
831 static 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 
863 std::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 
900 int 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 
907 template <int Dimensions>
908 static bool
909 Assemble(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 
1074  pmat = convert_to_hypre_matrix(
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 
1097 bool 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
void SetVerbose(const bool verbose)
Expose process details with verbose output.
int order_multiplier(const Geometry::Type type, const int dim)
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...
STL namespace.
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 AddMortarIntegrator(const std::shared_ptr< MortarIntegrator > &integrator)
This method must be called before Assemble or Transfer. It will assemble the operator in all intersec...
void source(const Vector &x, Vector &f)
Definition: ex25.cpp:620
void write(std::ostream &os, T value)
Write &#39;value&#39; to stream.
Definition: binaryio.hpp:30
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:912
ParMortarAssembler(const std::shared_ptr< ParFiniteElementSpace > &source, const std::shared_ptr< ParFiniteElementSpace > &destination)
constructs the object with source and destination spaces
T read(std::istream &is)
Read a value from the stream and return it.
Definition: binaryio.hpp:37
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3252
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)
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
std::shared_ptr< Cut > NewCut(const int dim)
Definition: cut.cpp:436
void SetMaxSolverIterations(const int max_solver_iterations)
Control the maximum numbers of conjugate gradients steps for mass matrix inversion.
string space
HYPRE_Int HYPRE_BigInt
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:116
double a
Definition: lissajous.cpp:41
int dim
Definition: ex24.cpp:53
bool Assemble(std::shared_ptr< HypreParMatrix > &B)
assembles the coupling matrix B. B : source -> destination If u is a coefficient associated with sour...
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
Class for parallel bilinear form.
void SetAssembleMassAndCouplingTogether(const bool value)
Control if the Mass matrix is computed together with the coupling operator every time.
Vector data type.
Definition: vector.hpp:60
Data type point element.
Definition: point.hpp:22
void pts(int iphi, int t, double x[])
void SetSolver(const std::shared_ptr< IterativeSolver > &solver)
Changes the solver to be used for solving the mass-matrix linear system.
RefCoord s[3]
Class for parallel grid function.
Definition: pgridfunc.hpp:32
bool Update()
assembles the various components necessary for the transfer. To be called before calling the Apply fu...