MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mortarassembler.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 "mortarassembler.hpp"
13 #include "../../general/tic_toc.hpp"
14 
15 #include "cut.hpp"
16 #include "transferutils.hpp"
17 
18 #include <cassert>
19 
20 // Moonolith includes
21 #include "moonolith_aabb.hpp"
22 #include "moonolith_serial_hash_grid.hpp"
23 #include "moonolith_stream_utils.hpp"
24 #include "par_moonolith_config.hpp"
25 
26 using namespace mfem::internal;
27 
28 namespace mfem
29 {
30 
31 struct MortarAssembler::Impl
32 {
33 public:
34  std::shared_ptr<FiniteElementSpace> source;
35  std::shared_ptr<FiniteElementSpace> destination;
36  std::vector<std::shared_ptr<MortarIntegrator>> integrators;
37  std::shared_ptr<SparseMatrix> coupling_matrix;
38  std::shared_ptr<SparseMatrix> mass_matrix;
39  bool verbose{false};
40  bool assemble_mass_and_coupling_together{true};
41  int max_solver_iterations{400};
42 
43  bool is_vector_fe() const
44  {
45  bool is_vector_fe = false;
46  for (auto i_ptr : integrators)
47  {
48  if (i_ptr->is_vector_fe())
49  {
50  is_vector_fe = true;
51  break;
52  }
53  }
54 
55  return is_vector_fe;
56  }
57 
58  BilinearFormIntegrator * new_mass_integrator() const
59  {
60  if (is_vector_fe())
61  {
62  return new VectorFEMassIntegrator();
63  }
64  else
65  {
66  return new MassIntegrator();
67  }
68  }
69 };
70 
72 
74 {
75  impl_->assemble_mass_and_coupling_together = value;
76 }
77 
78 void MortarAssembler::SetMaxSolverIterations(const int max_solver_iterations)
79 {
80  impl_->max_solver_iterations = max_solver_iterations;
81 }
82 
84  const std::shared_ptr<MortarIntegrator> &integrator)
85 {
86  impl_->integrators.push_back(integrator);
87 }
88 
89 void MortarAssembler::SetVerbose(const bool verbose)
90 {
91  impl_->verbose = verbose;
92 }
93 
94 template <int Dim>
95 void BuildBoxes(const Mesh &mesh,
96  std::vector<::moonolith::AABB<Dim, double>> &element_boxes)
97 {
98 #ifndef NDEBUG
99  const int dim = mesh.Dimension();
100  assert(dim == Dim);
101 #endif
102  element_boxes.resize(mesh.GetNE());
103 
105  for (int i = 0; i < mesh.GetNE(); ++i)
106  {
107  mesh.GetPointMatrix(i, pts);
108  MinCol(pts, &element_boxes[i].min_[0], false);
109  MaxCol(pts, &element_boxes[i].max_[0], false);
110  }
111 }
112 
113 bool HashGridDetectIntersections(const Mesh &src, const Mesh &dest,
114  std::vector<moonolith::Integer> &pairs)
115 {
116  const int dim = dest.Dimension();
117 
118  switch (dim)
119  {
120  case 1:
121  {
122  std::vector<::moonolith::AABB<1, double>> src_boxes, dest_boxes;
123  BuildBoxes(src, src_boxes);
124  BuildBoxes(dest, dest_boxes);
125 
126  ::moonolith::SerialHashGrid<1, double> grid;
127  return grid.detect(src_boxes, dest_boxes, pairs);
128  }
129  case 2:
130  {
131  std::vector<::moonolith::AABB<2, double>> src_boxes, dest_boxes;
132  BuildBoxes(src, src_boxes);
133  BuildBoxes(dest, dest_boxes);
134 
135  ::moonolith::SerialHashGrid<2, double> grid;
136  return grid.detect(src_boxes, dest_boxes, pairs);
137  }
138  case 3:
139  {
140  std::vector<::moonolith::AABB<3, double>> src_boxes, dest_boxes;
141  BuildBoxes(src, src_boxes);
142  BuildBoxes(dest, dest_boxes);
143 
144  ::moonolith::SerialHashGrid<3, double> grid;
145  return grid.detect(src_boxes, dest_boxes, pairs);
146  }
147  default:
148  {
149  assert(false);
150  return false;
151  }
152  }
153 }
154 
156  const std::shared_ptr<FiniteElementSpace> &source,
157  const std::shared_ptr<FiniteElementSpace> &destination)
158  : impl_(new Impl())
159 {
160  impl_->source = source;
161  impl_->destination = destination;
162 }
163 
164 int order_multiplier(const Geometry::Type type, const int dim)
165 {
166  return
167  (type == Geometry::TRIANGLE || type == Geometry::TETRAHEDRON ||
168  type == Geometry::SEGMENT)? 1 : dim;
169 }
170 
171 bool MortarAssembler::Assemble(std::shared_ptr<SparseMatrix> &B)
172 {
173  using namespace std;
174  const bool verbose = impl_->verbose;
175 
176  const auto &source_mesh = *impl_->source->GetMesh();
177  const auto &destination_mesh = *impl_->destination->GetMesh();
178 
179  int dim = source_mesh.Dimension();
180 
181  std::vector<::moonolith::Integer> pairs;
182  if (!HashGridDetectIntersections(source_mesh, destination_mesh, pairs))
183  {
184  return false;
185  }
186 
187  std::shared_ptr<Cut> cut = NewCut(dim);
188  if (!cut)
189  {
190  assert(false && "NOT Supported!");
191  return false;
192  }
193 
194 
195  IntegrationRule source_ir;
196  IntegrationRule destination_ir;
197 
198  int skip_zeros = 1;
199  B = make_shared<SparseMatrix>(impl_->destination->GetNDofs(),
200  impl_->source->GetNDofs());
201 
202 
203  std::unique_ptr<BilinearFormIntegrator> mass_integr(impl_->new_mass_integrator());
204 
205  if(impl_->assemble_mass_and_coupling_together) {
206  impl_->mass_matrix = make_shared<SparseMatrix>(impl_->destination->GetNDofs(), impl_->destination->GetNDofs());
207  }
208 
209 
210 
211  Array<int> source_vdofs, destination_vdofs;
212  DenseMatrix elemmat;
213  DenseMatrix cumulative_elemmat;
214  double local_element_matrices_sum = 0.0;
215 
216  long n_intersections = 0;
217  long n_candidates = 0;
218 
219  int max_q_order = 0;
220 
221  for (auto i_ptr : impl_->integrators)
222  {
223  max_q_order = std::max(i_ptr->GetQuadratureOrder(), max_q_order);
224  }
225 
226  bool intersected = false;
227  for (auto it = begin(pairs); it != end(pairs); /* inside */)
228  {
229  const int source_index = *it++;
230  const int destination_index = *it++;
231 
232  auto &source_fe = *impl_->source->GetFE(source_index);
233  auto &destination_fe = *impl_->destination->GetFE(destination_index);
234 
235  ElementTransformation &destination_Trans =
236  *impl_->destination->GetElementTransformation(destination_index);
237 
238  // Quadrature order mangling
239  int src_order_mult = order_multiplier(source_fe.GetGeomType(), dim);
240  int dest_order_mult = order_multiplier(destination_fe.GetGeomType(), dim);
241 
242  const int src_order = src_order_mult * source_fe.GetOrder();
243  const int dest_order = dest_order_mult * destination_fe.GetOrder();
244 
245  int contraction_order = src_order + dest_order;
246 
247  if(impl_->assemble_mass_and_coupling_together) {
248  contraction_order = std::max(contraction_order, 2 * dest_order);
249  }
250 
251  const int order = contraction_order + dest_order_mult * destination_Trans.OrderW() + max_q_order;
252 
253  // Update the quadrature rule in case it changed the order
254  cut->SetIntegrationOrder(order);
255 
256  n_candidates++;
257 
258  if (cut->BuildQuadrature(*impl_->source, source_index, *impl_->destination,
259  destination_index, source_ir, destination_ir))
260  {
261  impl_->source->GetElementVDofs(source_index, source_vdofs);
262  impl_->destination->GetElementVDofs(destination_index, destination_vdofs);
263 
264  ElementTransformation &source_Trans =
265  *impl_->source->GetElementTransformation(source_index);
266 
267  bool first = true;
268  for (auto i_ptr : impl_->integrators)
269  {
270  if (first)
271  {
272  i_ptr->AssembleElementMatrix(source_fe, source_ir, source_Trans,
273  destination_fe, destination_ir,
274  destination_Trans, cumulative_elemmat);
275  first = false;
276  }
277  else
278  {
279  i_ptr->AssembleElementMatrix(source_fe, source_ir, source_Trans,
280  destination_fe, destination_ir,
281  destination_Trans, elemmat);
282  cumulative_elemmat += elemmat;
283  }
284  }
285 
286  local_element_matrices_sum += Sum(cumulative_elemmat);
287 
288  B->AddSubMatrix(destination_vdofs, source_vdofs, cumulative_elemmat,
289  skip_zeros);
290 
291  if(impl_->assemble_mass_and_coupling_together) {
292  mass_integr->SetIntRule(&destination_ir);
293  mass_integr->AssembleElementMatrix(destination_fe, destination_Trans, elemmat);
294  impl_->mass_matrix->AddSubMatrix(destination_vdofs, destination_vdofs, elemmat, skip_zeros);
295  }
296 
297  intersected = true;
298  ++n_intersections;
299  }
300  }
301 
302  if (!intersected)
303  {
304  return false;
305  }
306 
307  B->Finalize();
308 
309  if(impl_->assemble_mass_and_coupling_together) {
310  impl_->mass_matrix->Finalize();
311  }
312 
313  if (verbose)
314  {
315  mfem::out << "local_element_matrices_sum: " << local_element_matrices_sum
316  << std::endl;
317  mfem::out << "B in R^(" << B->Height() << " x " << B->Width() << ")"
318  << std::endl;
319 
320  mfem::out << "n_intersections: " << n_intersections
321  << ", n_candidates: " << n_candidates << '\n';
322 
323  cut->Describe();
324  }
325 
326  return true;
327 }
328 
330  GridFunction &dest_fun)
331 {
332  return Update() && Apply(src_fun, dest_fun);
333 }
334 
336  GridFunction &dest_fun)
337 {
338  if (!impl_->coupling_matrix)
339  {
340  if (!Update())
341  {
342  return false;
343  }
344  }
345 
346  Vector temp(impl_->coupling_matrix->Height());
347  impl_->coupling_matrix->Mult(src_fun, temp);
348 
349  CGSolver Dinv;
350  Dinv.SetMaxIter(impl_->max_solver_iterations);
351 
352  if(impl_->verbose) {
353  Dinv.SetPrintLevel(3);
354  }
355 
356  Dinv.SetOperator(*impl_->mass_matrix);
357  Dinv.SetRelTol(1e-6);
358  Dinv.SetMaxIter(80);
359  Dinv.Mult(temp, dest_fun);
360  return true;
361 }
362 
364 {
365  using namespace std;
366  const bool verbose = impl_->verbose;
367 
368  StopWatch chrono;
369 
370  if (verbose)
371  {
372  mfem::out << "\nAssembling coupling operator..." << endl;
373  }
374 
375  chrono.Start();
376 
377  if (!Assemble(impl_->coupling_matrix))
378  {
379  return false;
380  }
381 
382  chrono.Stop();
383  if (verbose)
384  {
385  mfem::out << "Done. time: ";
386  mfem::out << chrono.RealTime() << " seconds" << endl;
387  }
388 
389  if(!impl_->assemble_mass_and_coupling_together) {
390  BilinearForm b_form(impl_->destination.get());
391 
392  b_form.AddDomainIntegrator(impl_->new_mass_integrator());
393 
394  b_form.Assemble();
395  b_form.Finalize();
396 
397  impl_->mass_matrix = std::shared_ptr<SparseMatrix>(b_form.LoseMat());
398  }
399 
400  if (verbose)
401  {
402  Vector brs(impl_->coupling_matrix->Height());
403  impl_->coupling_matrix->GetRowSums(brs);
404 
405  Vector drs(impl_->mass_matrix->Height());
406  impl_->mass_matrix->GetRowSums(drs);
407 
408  mfem::out << "sum(B): " << brs.Sum() << std::endl;
409  mfem::out << "sum(D): " << drs.Sum() << std::endl;
410  }
411 
412  return true;
413 }
414 
415 } // namespace mfem
void GetPointMatrix(int i, DenseMatrix &pointmat) const
Definition: mesh.cpp:6260
Conjugate gradient method.
Definition: solvers.hpp:465
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int order_multiplier(const Geometry::Type type, const int dim)
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
bool Apply(const GridFunction &src_fun, GridFunction &dest_fun)
transfer a function from source to destination. It requires that the Update function is called before...
void SetVerbose(const bool verbose)
Expose process details with verbose output.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
bool HashGridDetectIntersections(const Mesh &src, const Mesh &dest, std::vector< moonolith::Integer > &pairs)
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
bool Update()
assembles the various components necessary for the transfer. To be called before calling the Apply fu...
void source(const Vector &x, Vector &f)
Definition: ex25.cpp:620
MortarAssembler(const std::shared_ptr< FiniteElementSpace > &source, const std::shared_ptr< FiniteElementSpace > &destination)
constructs the object with source and destination spaces
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
bool Transfer(const GridFunction &src_fun, GridFunction &dest_fun)
transfer a function from source to destination. if the transfer is to be performed multiple times use...
Timing object.
Definition: tic_toc.hpp:34
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
int Dimension() const
Definition: mesh.hpp:1006
void SetAssembleMassAndCouplingTogether(const bool value)
Control if the Mass matrix is computed together with the coupling operator every time.
void SetRelTol(double rtol)
Definition: solvers.hpp:198
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.
void BuildBoxes(const Mesh &mesh, std::vector<::moonolith::AABB< Dim, double >> &element_boxes)
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
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...
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
bool Assemble(std::shared_ptr< SparseMatrix > &B)
assembles the coupling matrix B. B : source -&gt; destination If u is a coefficient associated with sour...
Vector data type.
Definition: vector.hpp:60
void pts(int iphi, int t, double x[])
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