Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #ifdef MOOSE_MFEM_ENABLED
11 :
12 : #include "PatchedMixedBilinearForm.h"
13 :
14 : namespace Moose::MFEM
15 : {
16 :
17 : void
18 173 : ParMixedBilinearForm::Assemble(int skip_zeros)
19 : {
20 173 : if (interior_face_integs.Size())
21 : {
22 0 : trial_pfes->ExchangeFaceNbrData();
23 0 : test_pfes->ExchangeFaceNbrData();
24 0 : if (!ext && mat == NULL)
25 : {
26 0 : pAllocMat();
27 : }
28 : }
29 :
30 : // Call SubMeshTolerantAssemble if trial_pfes is defined on a submesh of the mesh
31 : // on which test_pfes is defined, else default to MixedBilinearForm::Assemble
32 : // as usual
33 184 : if ((trial_pfes->GetParMesh() != test_pfes->GetParMesh()) &&
34 11 : mfem::ParSubMesh::IsParSubMesh(trial_pfes->GetParMesh(), test_pfes->GetParMesh()))
35 11 : SubMeshTolerantAssemble(skip_zeros);
36 : else
37 162 : MixedBilinearForm::Assemble(skip_zeros);
38 :
39 173 : if (!ext && interior_face_integs.Size() > 0)
40 : {
41 0 : AssembleSharedFaces(skip_zeros);
42 : }
43 173 : }
44 :
45 : void
46 11 : ParMixedBilinearForm::SubMeshTolerantAssemble(int skip_zeros)
47 : {
48 11 : if (ext)
49 : {
50 0 : ext->Assemble();
51 0 : return;
52 : }
53 :
54 : mfem::ElementTransformation * eltrans;
55 11 : mfem::DenseMatrix elmat;
56 :
57 : // For usual mfem::MixedBilinearFormAssemble, the iterated FESpace is test_fes
58 11 : mfem::FiniteElementSpace & iterated_fes = *trial_fes;
59 11 : mfem::Mesh * mesh = iterated_fes.GetMesh();
60 11 : const auto * trial_submesh = static_cast<const mfem::ParSubMesh *>(trial_pfes->GetParMesh());
61 11 : const auto * trial_parent_mesh = trial_submesh->GetParent();
62 11 : const auto & test_element_ids = trial_submesh->GetParentElementIDMap();
63 11 : const auto & test_face_ids = trial_submesh->GetParentFaceIDMap();
64 11 : const auto & test_edge_ids = trial_submesh->GetParentEdgeIDMap();
65 11 : const auto & test_vertex_ids = trial_submesh->GetParentVertexIDMap();
66 11 : const auto & test_face_to_be = trial_parent_mesh->GetFaceToBdrElMap();
67 11 : const int num_elem = iterated_fes.GetNE();
68 11 : const int num_boundary_elem = iterated_fes.GetNBE();
69 :
70 11 : if (mat == NULL)
71 : {
72 11 : mat = new mfem::SparseMatrix(height, width);
73 : }
74 :
75 11 : if (domain_integs.Size())
76 : {
77 22 : for (int k = 0; k < domain_integs.Size(); k++)
78 : {
79 11 : if (domain_integs_marker[k] != NULL)
80 : {
81 0 : MFEM_VERIFY(domain_integs_marker[k]->Size() ==
82 : (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
83 : "invalid element marker for domain integrator #" << k
84 : << ", counting from zero");
85 : }
86 : }
87 :
88 11 : mfem::DofTransformation dom_dof_trans, ran_dof_trans;
89 1019 : for (int i = 0; i < num_elem; i++)
90 : {
91 1008 : const int test_elem_id = test_element_ids[i];
92 1008 : const int elem_attr = mesh->GetAttribute(i);
93 1008 : trial_fes->GetElementVDofs(i, trial_vdofs, dom_dof_trans);
94 1008 : test_fes->GetElementVDofs(test_elem_id, test_vdofs, ran_dof_trans);
95 1008 : eltrans = iterated_fes.GetElementTransformation(i);
96 :
97 1008 : elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
98 1008 : elmat = 0.0;
99 2016 : for (int k = 0; k < domain_integs.Size(); k++)
100 : {
101 1008 : if (domain_integs_marker[k] == NULL || (*(domain_integs_marker[k]))[elem_attr - 1] == 1)
102 : {
103 1008 : domain_integs[k]->AssembleElementMatrix2(
104 1008 : *trial_fes->GetFE(i), *test_fes->GetFE(test_elem_id), *eltrans, elemmat);
105 1008 : elmat += elemmat;
106 : }
107 : }
108 1008 : TransformDual(ran_dof_trans, dom_dof_trans, elmat);
109 1008 : mat->AddSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
110 : }
111 11 : }
112 :
113 11 : if (boundary_integs.Size())
114 : {
115 : // Which boundary attributes need to be processed?
116 0 : mfem::Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ? mesh->bdr_attributes.Max() : 0);
117 0 : bdr_attr_marker = 0;
118 0 : for (int k = 0; k < boundary_integs.Size(); k++)
119 : {
120 0 : if (boundary_integs_marker[k] == NULL)
121 : {
122 0 : bdr_attr_marker = 1;
123 0 : break;
124 : }
125 0 : mfem::Array<int> & bdr_marker = *boundary_integs_marker[k];
126 : MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
127 : "invalid boundary marker for boundary integrator #" << k
128 : << ", counting from zero");
129 0 : for (int i = 0; i < bdr_attr_marker.Size(); i++)
130 : {
131 0 : bdr_attr_marker[i] |= bdr_marker[i];
132 : }
133 : }
134 :
135 0 : mfem::DofTransformation dom_dof_trans, ran_dof_trans;
136 0 : for (int i = 0; i < num_boundary_elem; i++)
137 : {
138 0 : const int iface = mesh->GetBdrElementFaceIndex(i);
139 0 : const int test_face_id = mesh->Dimension() == 3 ? test_face_ids[iface]
140 0 : : mesh->Dimension() == 2 ? test_edge_ids[iface]
141 0 : : test_vertex_ids[iface];
142 0 : const int test_bdr_elem_id = test_face_to_be[test_face_id];
143 0 : const int bdr_attr = mesh->GetBdrAttribute(i);
144 0 : if (bdr_attr_marker[bdr_attr - 1] == 0)
145 : {
146 0 : continue;
147 : }
148 0 : MFEM_VERIFY(test_bdr_elem_id != -1,
149 : "Cannot assemble a mixed boundary integrator from submesh boundary element "
150 : << i << " because it does not correspond to a parent boundary element.");
151 :
152 0 : trial_fes->GetBdrElementVDofs(i, trial_vdofs, dom_dof_trans);
153 0 : test_fes->GetBdrElementVDofs(test_bdr_elem_id, test_vdofs, ran_dof_trans);
154 0 : eltrans = iterated_fes.GetBdrElementTransformation(i);
155 :
156 0 : elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
157 0 : elmat = 0.0;
158 0 : for (int k = 0; k < boundary_integs.Size(); k++)
159 : {
160 0 : if (boundary_integs_marker[k] && (*boundary_integs_marker[k])[bdr_attr - 1] == 0)
161 : {
162 0 : continue;
163 : }
164 :
165 0 : boundary_integs[k]->AssembleElementMatrix2(
166 0 : *trial_fes->GetBE(i), *test_fes->GetBE(test_bdr_elem_id), *eltrans, elemmat);
167 0 : elmat += elemmat;
168 : }
169 0 : TransformDual(ran_dof_trans, dom_dof_trans, elmat);
170 0 : mat->AddSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
171 : }
172 0 : }
173 :
174 11 : if (interior_face_integs.Size())
175 : {
176 : mfem::FaceElementTransformations * ftr;
177 0 : mfem::Array<int> trial_vdofs2, test_vdofs2;
178 : const mfem::FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
179 :
180 0 : int nfaces = mesh->GetNumFaces();
181 0 : for (int i = 0; i < nfaces; i++)
182 : {
183 0 : ftr = mesh->GetInteriorFaceTransformations(i);
184 0 : if (ftr != NULL)
185 : {
186 0 : const int test_elem_1 = test_element_ids[ftr->Elem1No];
187 0 : trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
188 0 : test_fes->GetElementVDofs(test_elem_1, test_vdofs);
189 0 : trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
190 0 : test_fe1 = test_fes->GetFE(test_elem_1);
191 0 : if (ftr->Elem2No >= 0)
192 : {
193 0 : const int test_elem_2 = test_element_ids[ftr->Elem2No];
194 0 : trial_fes->GetElementVDofs(ftr->Elem2No, trial_vdofs2);
195 0 : test_fes->GetElementVDofs(test_elem_2, test_vdofs2);
196 0 : trial_vdofs.Append(trial_vdofs2);
197 0 : test_vdofs.Append(test_vdofs2);
198 0 : trial_fe2 = trial_fes->GetFE(ftr->Elem2No);
199 0 : test_fe2 = test_fes->GetFE(test_elem_2);
200 : }
201 : else
202 : {
203 : // The test_fe2 object is really a dummy and not used on the
204 : // boundaries, but we can't dereference a NULL pointer, and we don't
205 : // want to actually make a fake element.
206 0 : trial_fe2 = trial_fe1;
207 0 : test_fe2 = test_fe1;
208 : }
209 0 : for (int k = 0; k < interior_face_integs.Size(); k++)
210 : {
211 0 : interior_face_integs[k]->AssembleFaceMatrix(
212 0 : *trial_fe1, *test_fe1, *trial_fe2, *test_fe2, *ftr, elemmat);
213 0 : mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
214 : }
215 : }
216 : }
217 0 : }
218 :
219 11 : if (boundary_face_integs.Size())
220 : {
221 : mfem::FaceElementTransformations * ftr;
222 0 : mfem::Array<int> tr_vdofs2, te_vdofs2;
223 : const mfem::FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
224 :
225 : // Which boundary attributes need to be processed?
226 0 : mfem::Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ? mesh->bdr_attributes.Max() : 0);
227 0 : bdr_attr_marker = 0;
228 0 : for (int k = 0; k < boundary_face_integs.Size(); k++)
229 : {
230 0 : if (boundary_face_integs_marker[k] == NULL)
231 : {
232 0 : bdr_attr_marker = 1;
233 0 : break;
234 : }
235 0 : mfem::Array<int> & bdr_marker = *boundary_face_integs_marker[k];
236 : MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
237 : "invalid boundary marker for boundary face integrator #"
238 : << k << ", counting from zero");
239 0 : for (int i = 0; i < bdr_attr_marker.Size(); i++)
240 : {
241 0 : bdr_attr_marker[i] |= bdr_marker[i];
242 : }
243 : }
244 :
245 0 : mfem::DofTransformation dom_dof_trans, ran_dof_trans;
246 0 : for (int i = 0; i < trial_fes->GetNBE(); i++)
247 : {
248 0 : const int bdr_attr = mesh->GetBdrAttribute(i);
249 0 : if (bdr_attr_marker[bdr_attr - 1] == 0)
250 : {
251 0 : continue;
252 : }
253 :
254 0 : ftr = mesh->GetBdrFaceTransformations(i);
255 0 : if (ftr != NULL)
256 : {
257 0 : const int test_elem_1 = test_element_ids[ftr->Elem1No];
258 0 : trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs, dom_dof_trans);
259 0 : test_fes->GetElementVDofs(test_elem_1, test_vdofs, ran_dof_trans);
260 0 : trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
261 0 : test_fe1 = test_fes->GetFE(test_elem_1);
262 : // The test_fe2 object is really a dummy and not used on the
263 : // boundaries, but we can't dereference a NULL pointer, and we don't
264 : // want to actually make a fake element.
265 0 : trial_fe2 = trial_fe1;
266 0 : test_fe2 = test_fe1;
267 0 : for (int k = 0; k < boundary_face_integs.Size(); k++)
268 : {
269 0 : if (boundary_face_integs_marker[k] &&
270 0 : (*boundary_face_integs_marker[k])[bdr_attr - 1] == 0)
271 : {
272 0 : continue;
273 : }
274 :
275 0 : boundary_face_integs[k]->AssembleFaceMatrix(
276 0 : *trial_fe1, *test_fe1, *trial_fe2, *test_fe2, *ftr, elemmat);
277 0 : TransformDual(ran_dof_trans, dom_dof_trans, elemmat);
278 0 : mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
279 : }
280 : }
281 : }
282 0 : }
283 :
284 11 : if (trace_face_integs.Size())
285 : {
286 : mfem::FaceElementTransformations * ftr;
287 0 : mfem::Array<int> test_vdofs2;
288 : const mfem::FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
289 :
290 0 : int nfaces = mesh->GetNumFaces();
291 0 : for (int i = 0; i < nfaces; i++)
292 : {
293 0 : ftr = mesh->GetFaceElementTransformations(i);
294 0 : const int test_elem_1 = test_element_ids[ftr->Elem1No];
295 0 : trial_fes->GetFaceVDofs(i, trial_vdofs);
296 0 : test_fes->GetElementVDofs(test_elem_1, test_vdofs);
297 0 : trial_face_fe = trial_fes->GetFaceElement(i);
298 0 : test_fe1 = test_fes->GetFE(test_elem_1);
299 0 : if (ftr->Elem2No >= 0)
300 : {
301 0 : const int test_elem_2 = test_element_ids[ftr->Elem2No];
302 0 : test_fes->GetElementVDofs(test_elem_2, test_vdofs2);
303 0 : test_vdofs.Append(test_vdofs2);
304 0 : test_fe2 = test_fes->GetFE(test_elem_2);
305 : }
306 : else
307 : {
308 : // The test_fe2 object is really a dummy and not used on the
309 : // boundaries, but we can't dereference a NULL pointer, and we don't
310 : // want to actually make a fake element.
311 0 : test_fe2 = test_fe1;
312 : }
313 0 : for (int k = 0; k < trace_face_integs.Size(); k++)
314 : {
315 0 : trace_face_integs[k]->AssembleFaceMatrix(
316 0 : *trial_face_fe, *test_fe1, *test_fe2, *ftr, elemmat);
317 0 : mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
318 : }
319 : }
320 0 : }
321 :
322 11 : if (boundary_trace_face_integs.Size())
323 : {
324 : mfem::FaceElementTransformations * ftr;
325 0 : mfem::Array<int> te_vdofs2;
326 : const mfem::FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
327 :
328 : // Which boundary attributes need to be processed?
329 0 : mfem::Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ? mesh->bdr_attributes.Max() : 0);
330 0 : bdr_attr_marker = 0;
331 0 : for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
332 : {
333 0 : if (boundary_trace_face_integs_marker[k] == NULL)
334 : {
335 0 : bdr_attr_marker = 1;
336 0 : break;
337 : }
338 0 : mfem::Array<int> & bdr_marker = *boundary_trace_face_integs_marker[k];
339 : MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
340 : "invalid boundary marker for boundary trace face"
341 : "integrator #"
342 : << k << ", counting from zero");
343 0 : for (int i = 0; i < bdr_attr_marker.Size(); i++)
344 : {
345 0 : bdr_attr_marker[i] |= bdr_marker[i];
346 : }
347 : }
348 :
349 0 : for (int i = 0; i < trial_fes->GetNBE(); i++)
350 : {
351 0 : const int bdr_attr = mesh->GetBdrAttribute(i);
352 0 : if (bdr_attr_marker[bdr_attr - 1] == 0)
353 : {
354 0 : continue;
355 : }
356 :
357 0 : ftr = mesh->GetBdrFaceTransformations(i);
358 0 : if (ftr)
359 : {
360 0 : const int iface = mesh->GetBdrElementFaceIndex(i);
361 0 : const int test_elem_1 = test_element_ids[ftr->Elem1No];
362 0 : trial_fes->GetFaceVDofs(iface, trial_vdofs);
363 0 : test_fes->GetElementVDofs(test_elem_1, test_vdofs);
364 0 : trial_face_fe = trial_fes->GetFaceElement(iface);
365 0 : test_fe1 = test_fes->GetFE(test_elem_1);
366 : // The test_fe2 object is really a dummy and not used on the
367 : // boundaries, but we can't dereference a NULL pointer, and we don't
368 : // want to actually make a fake element.
369 0 : test_fe2 = test_fe1;
370 0 : for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
371 : {
372 0 : if (boundary_trace_face_integs_marker[k] &&
373 0 : (*boundary_trace_face_integs_marker[k])[bdr_attr - 1] == 0)
374 : {
375 0 : continue;
376 : }
377 :
378 0 : boundary_trace_face_integs[k]->AssembleFaceMatrix(
379 0 : *trial_face_fe, *test_fe1, *test_fe2, *ftr, elemmat);
380 0 : mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
381 : }
382 : }
383 : }
384 0 : }
385 11 : }
386 :
387 : } // namespace Moose::MFEM
388 :
389 : #endif
|
