MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mortarassembler.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
12#include "mortarassembler.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
26using namespace mfem::internal;
27
28namespace mfem
29{
30
31struct MortarAssembler::Impl
32{
33public:
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
78void 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
89void MortarAssembler::SetVerbose(const bool verbose)
90{
91 impl_->verbose = verbose;
92}
93
94template <int Dim>
95void 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
113bool 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
164int 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
171bool 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
Abstract base class BilinearFormIntegrator.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
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.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
SparseMatrix * LoseMat()
Nullifies the internal matrix and returns a pointer to it. Used for transferring ownership.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Mesh data type.
Definition mesh.hpp:56
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetPointMatrix(int i, DenseMatrix &pointmat) const
Definition mesh.cpp:7326
void SetAssembleMassAndCouplingTogether(const bool value)
Control if the Mass matrix is computed together with the coupling operator every time.
void SetMaxSolverIterations(const int max_solver_iterations)
Control the maximum numbers of conjugate gradients steps for mass matrix inversion.
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...
MortarAssembler(const std::shared_ptr< FiniteElementSpace > &source, const std::shared_ptr< FiniteElementSpace > &destination)
constructs the object with source and destination spaces
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
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...
bool Assemble(std::shared_ptr< SparseMatrix > &B)
assembles the coupling matrix B. B : source -> destination If u is a coefficient associated with sour...
bool Update()
assembles the various components necessary for the transfer. To be called before calling the Apply fu...
Timing object.
Definition tic_toc.hpp:36
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Stop()
Stop the stopwatch.
Definition tic_toc.cpp:422
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
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 BuildBoxes(const Mesh &mesh, std::vector<::moonolith::AABB< Dim, double > > &element_boxes)
bool HashGridDetectIntersections(const Mesh &src, const Mesh &dest, std::vector< moonolith::Integer > &pairs)
int order_multiplier(const Geometry::Type type, const int dim)
void pts(int iphi, int t, real_t x[])