Skip to content

Commit 8bd7b33

Browse files
nmnobrenmnobre
committed
Use Lagrange mult to constrain the scalar field in the div-grad example
1 parent 788b669 commit 8bd7b33

2 files changed

Lines changed: 94 additions & 81 deletions

File tree

‎doc/html/examples/vector_fe_ex6.html‎

Lines changed: 72 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -69,30 +69,30 @@
6969
n_systems()=1
7070
System #0, "DivGrad"
7171
Type "LinearImplicit"
72-
Variables="u" "p"
73-
Finite Element Types="RAVIART_THOMAS" "MONOMIAL"
74-
Approximation Orders="FIRST" "CONSTANT"
75-
n_dofs()=1155
76-
n_local_dofs()=1155
77-
max(n_local_dofs())=1155
72+
Variables="u" "p" "l"
73+
Finite Element Types="RAVIART_THOMAS" "MONOMIAL" "SCALAR"
74+
Approximation Orders="FIRST" "CONSTANT" "FIRST"
75+
n_dofs()=1156
76+
n_local_dofs()=1156
77+
max(n_local_dofs())=1156
7878
n_constrained_dofs()=0
7979
n_local_constrained_dofs()=0
80-
max(local unconstrained dofs)=1155
80+
max(local unconstrained dofs)=1156
8181
n_vectors()=1
8282
n_matrices()=1
8383
DofMap Sparsity
84-
Average On-Processor Bandwidth <= 5.67532
84+
Average On-Processor Bandwidth <= 7.66955
8585
Average Off-Processor Bandwidth <= 0
86-
Maximum On-Processor Bandwidth <= 7
86+
Maximum On-Processor Bandwidth <= 1156
8787
Maximum Off-Processor Bandwidth <= 0
8888
DofMap Constraints
8989
Number of DoF Constraints = 0
9090

9191
WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Sep 11 2023 at 16:24:29 ***
9292
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Sep 11 2023 at 16:24:29 ***
93-
L2 error is: 0.134294
94-
HDiv semi-norm error is: 0.344131
95-
HDiv error is: 0.369406
93+
L2 error is: 0.13433
94+
HDiv semi-norm error is: 0.344125
95+
HDiv error is: 0.369414
9696

9797
***************************************************************
9898
* Done Running Example vector_fe_ex6:
@@ -124,30 +124,30 @@
124124
n_systems()=1
125125
System #0, "DivGrad"
126126
Type "LinearImplicit"
127-
Variables="u" "p"
128-
Finite Element Types="RAVIART_THOMAS" "MONOMIAL"
129-
Approximation Orders="FIRST" "CONSTANT"
130-
n_dofs()=1155
131-
n_local_dofs()=1155
132-
max(n_local_dofs())=1155
127+
Variables="u" "p" "l"
128+
Finite Element Types="RAVIART_THOMAS" "MONOMIAL" "SCALAR"
129+
Approximation Orders="FIRST" "CONSTANT" "FIRST"
130+
n_dofs()=1156
131+
n_local_dofs()=1156
132+
max(n_local_dofs())=1156
133133
n_constrained_dofs()=0
134134
n_local_constrained_dofs()=0
135-
max(local unconstrained dofs)=1155
135+
max(local unconstrained dofs)=1156
136136
n_vectors()=1
137137
n_matrices()=1
138138
DofMap Sparsity
139-
Average On-Processor Bandwidth <= 5.67532
139+
Average On-Processor Bandwidth <= 7.66955
140140
Average Off-Processor Bandwidth <= 0
141-
Maximum On-Processor Bandwidth <= 7
141+
Maximum On-Processor Bandwidth <= 1156
142142
Maximum Off-Processor Bandwidth <= 0
143143
DofMap Constraints
144144
Number of DoF Constraints = 0
145145

146146
WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Sep 11 2023 at 16:24:29 ***
147147
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Sep 11 2023 at 16:24:29 ***
148-
L2 error is: 0.134294
149-
HDiv semi-norm error is: 0.344131
150-
HDiv error is: 0.369406
148+
L2 error is: 0.13433
149+
HDiv semi-norm error is: 0.344125
150+
HDiv error is: 0.369414
151151

152152
***************************************************************
153153
* Done Running Example vector_fe_ex6:
@@ -179,30 +179,30 @@
179179
n_systems()=1
180180
System #0, "DivGrad"
181181
Type "LinearImplicit"
182-
Variables="u" "p"
183-
Finite Element Types="RAVIART_THOMAS" "MONOMIAL"
184-
Approximation Orders="FIRST" "CONSTANT"
185-
n_dofs()=705
186-
n_local_dofs()=705
187-
max(n_local_dofs())=705
182+
Variables="u" "p" "l"
183+
Finite Element Types="RAVIART_THOMAS" "MONOMIAL" "SCALAR"
184+
Approximation Orders="FIRST" "CONSTANT" "FIRST"
185+
n_dofs()=706
186+
n_local_dofs()=706
187+
max(n_local_dofs())=706
188188
n_constrained_dofs()=0
189189
n_local_constrained_dofs()=0
190-
max(local unconstrained dofs)=705
190+
max(local unconstrained dofs)=706
191191
n_vectors()=1
192192
n_matrices()=1
193193
DofMap Sparsity
194-
Average On-Processor Bandwidth <= 7.38298
194+
Average On-Processor Bandwidth <= 9.3711
195195
Average Off-Processor Bandwidth <= 0
196-
Maximum On-Processor Bandwidth <= 9
196+
Maximum On-Processor Bandwidth <= 706
197197
Maximum Off-Processor Bandwidth <= 0
198198
DofMap Constraints
199199
Number of DoF Constraints = 0
200200

201201
WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Sep 11 2023 at 16:24:29 ***
202202
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Sep 11 2023 at 16:24:29 ***
203-
L2 error is: 0.134405
204-
HDiv semi-norm error is: 0.421309
205-
HDiv error is: 0.442228
203+
L2 error is: 0.13447
204+
HDiv semi-norm error is: 0.4213
205+
HDiv error is: 0.44224
206206

207207
***************************************************************
208208
* Done Running Example vector_fe_ex6:
@@ -234,30 +234,30 @@
234234
n_systems()=1
235235
System #0, "DivGrad"
236236
Type "LinearImplicit"
237-
Variables="u" "p"
238-
Finite Element Types="RAVIART_THOMAS" "MONOMIAL"
239-
Approximation Orders="FIRST" "CONSTANT"
240-
n_dofs()=705
241-
n_local_dofs()=705
242-
max(n_local_dofs())=705
237+
Variables="u" "p" "l"
238+
Finite Element Types="RAVIART_THOMAS" "MONOMIAL" "SCALAR"
239+
Approximation Orders="FIRST" "CONSTANT" "FIRST"
240+
n_dofs()=706
241+
n_local_dofs()=706
242+
max(n_local_dofs())=706
243243
n_constrained_dofs()=0
244244
n_local_constrained_dofs()=0
245-
max(local unconstrained dofs)=705
245+
max(local unconstrained dofs)=706
246246
n_vectors()=1
247247
n_matrices()=1
248248
DofMap Sparsity
249-
Average On-Processor Bandwidth <= 7.38298
249+
Average On-Processor Bandwidth <= 9.3711
250250
Average Off-Processor Bandwidth <= 0
251-
Maximum On-Processor Bandwidth <= 9
251+
Maximum On-Processor Bandwidth <= 706
252252
Maximum Off-Processor Bandwidth <= 0
253253
DofMap Constraints
254254
Number of DoF Constraints = 0
255255

256256
WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Sep 11 2023 at 16:24:29 ***
257257
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Sep 11 2023 at 16:24:29 ***
258-
L2 error is: 0.134405
259-
HDiv semi-norm error is: 0.421309
260-
HDiv error is: 0.442228
258+
L2 error is: 0.13447
259+
HDiv semi-norm error is: 0.4213
260+
HDiv error is: 0.44224
261261

262262
***************************************************************
263263
* Done Running Example vector_fe_ex6:
@@ -289,30 +289,30 @@
289289
n_systems()=1
290290
System #0, "DivGrad"
291291
Type "LinearImplicit"
292-
Variables="u" "p"
293-
Finite Element Types="RAVIART_THOMAS" "MONOMIAL"
294-
Approximation Orders="FIRST" "CONSTANT"
295-
n_dofs()=15984
296-
n_local_dofs()=15984
297-
max(n_local_dofs())=15984
292+
Variables="u" "p" "l"
293+
Finite Element Types="RAVIART_THOMAS" "MONOMIAL" "SCALAR"
294+
Approximation Orders="FIRST" "CONSTANT" "FIRST"
295+
n_dofs()=15985
296+
n_local_dofs()=15985
297+
max(n_local_dofs())=15985
298298
n_constrained_dofs()=0
299299
n_local_constrained_dofs()=0
300-
max(local unconstrained dofs)=15984
300+
max(local unconstrained dofs)=15985
301301
n_vectors()=1
302302
n_matrices()=1
303303
DofMap Sparsity
304-
Average On-Processor Bandwidth <= 7.48649
304+
Average On-Processor Bandwidth <= 9.48596
305305
Average Off-Processor Bandwidth <= 0
306-
Maximum On-Processor Bandwidth <= 9
306+
Maximum On-Processor Bandwidth <= 15985
307307
Maximum Off-Processor Bandwidth <= 0
308308
DofMap Constraints
309309
Number of DoF Constraints = 0
310310

311311
WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Sep 11 2023 at 16:24:29 ***
312312
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Sep 11 2023 at 16:24:29 ***
313-
L2 error is: 0.33466
314-
HDiv semi-norm error is: 0.837163
315-
HDiv error is: 0.901576
313+
L2 error is: 0.334676
314+
HDiv semi-norm error is: 0.83716
315+
HDiv error is: 0.90158
316316

317317
***************************************************************
318318
* Done Running Example vector_fe_ex6:
@@ -344,30 +344,30 @@
344344
n_systems()=1
345345
System #0, "DivGrad"
346346
Type "LinearImplicit"
347-
Variables="u" "p"
348-
Finite Element Types="RAVIART_THOMAS" "MONOMIAL"
349-
Approximation Orders="FIRST" "CONSTANT"
350-
n_dofs()=14175
351-
n_local_dofs()=14175
352-
max(n_local_dofs())=14175
347+
Variables="u" "p" "l"
348+
Finite Element Types="RAVIART_THOMAS" "MONOMIAL" "SCALAR"
349+
Approximation Orders="FIRST" "CONSTANT" "FIRST"
350+
n_dofs()=14176
351+
n_local_dofs()=14176
352+
max(n_local_dofs())=14176
353353
n_constrained_dofs()=0
354354
n_local_constrained_dofs()=0
355-
max(local unconstrained dofs)=14175
355+
max(local unconstrained dofs)=14176
356356
n_vectors()=1
357357
n_matrices()=1
358358
DofMap Sparsity
359-
Average On-Processor Bandwidth <= 11
359+
Average On-Processor Bandwidth <= 12.9992
360360
Average Off-Processor Bandwidth <= 0
361-
Maximum On-Processor Bandwidth <= 13
361+
Maximum On-Processor Bandwidth <= 14176
362362
Maximum Off-Processor Bandwidth <= 0
363363
DofMap Constraints
364364
Number of DoF Constraints = 0
365365

366366
WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Sep 11 2023 at 16:24:29 ***
367367
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Sep 11 2023 at 16:24:29 ***
368-
L2 error is: 0.232445
369-
HDiv semi-norm error is: 0.773276
370-
HDiv error is: 0.807457
368+
L2 error is: 0.232486
369+
HDiv semi-norm error is: 0.773271
370+
HDiv error is: 0.807464
371371

372372
***************************************************************
373373
* Done Running Example vector_fe_ex6:

‎examples/vector_fe/vector_fe_ex6/vector_fe_ex6.C‎

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -145,6 +145,9 @@ int main (int argc, char ** argv)
145145
system.add_variable("u", FIRST, RAVIART_THOMAS);
146146
system.add_variable("p", CONSTANT, MONOMIAL);
147147

148+
// Add a scalar Lagrange multiplier to enforce our constraint
149+
system.add_variable("l", FIRST, SCALAR);
150+
148151
// Give the system a pointer to the matrix assembly
149152
// function. This will be called when needed by the library.
150153
system.attach_assemble_function(assemble_divgrad);
@@ -310,6 +313,7 @@ void assemble_divgrad(EquationSystems & es,
310313
std::vector<dof_id_type> dof_indices;
311314
std::vector<dof_id_type> vector_dof_indices;
312315
std::vector<dof_id_type> scalar_dof_indices;
316+
std::vector<dof_id_type> lambda_dof_indices;
313317

314318
// The global system matrix
315319
SparseMatrix<Number> & matrix = system.get_system_matrix();
@@ -336,6 +340,7 @@ void assemble_divgrad(EquationSystems & es,
336340
dof_map.dof_indices (elem, dof_indices);
337341
dof_map.dof_indices (elem, vector_dof_indices, system.variable_number("u"));
338342
dof_map.dof_indices (elem, scalar_dof_indices, system.variable_number("p"));
343+
dof_map.dof_indices (elem, lambda_dof_indices, system.variable_number("l"));
339344

340345
// Cache the number of degrees of freedom, in total and for each
341346
// variable, on this element, for use as array and loop bounds later.
@@ -348,6 +353,8 @@ void assemble_divgrad(EquationSystems & es,
348353
cast_int<unsigned int>(vector_dof_indices.size());
349354
const unsigned int scalar_n_dofs =
350355
cast_int<unsigned int>(scalar_dof_indices.size());
356+
const unsigned int lambda_n_dofs =
357+
cast_int<unsigned int>(lambda_dof_indices.size());
351358

352359
// Compute the element-specific data for the current
353360
// element. This involves computing the location of the
@@ -359,7 +366,7 @@ void assemble_divgrad(EquationSystems & es,
359366
// The total number of degrees of freedom is just the sum of the number
360367
// of degrees of freedom per variable. We should also have the same
361368
// number of degrees of freedom as shape functions for each variable.
362-
libmesh_assert_equal_to (n_dofs, vector_n_dofs + scalar_n_dofs);
369+
libmesh_assert_equal_to (n_dofs, vector_n_dofs + scalar_n_dofs + lambda_n_dofs);
363370
libmesh_assert_equal_to (vector_n_dofs, vector_phi.size());
364371
libmesh_assert_equal_to (scalar_n_dofs, scalar_phi.size());
365372

@@ -452,20 +459,26 @@ void assemble_divgrad(EquationSystems & es,
452459
else if (dim == 3)
453460
scalar_value = DivGradExactSolution().scalar(x, y, z);
454461

455-
// The lower-right block involves a double loop to integrate the
456-
// scalar test functions (k) against the
457-
// scalar trial functions (l).
462+
// A double loop to integrate the
463+
// scalar test functions (k) against the Lagrange dof (n).
458464
for (unsigned int k = 0; k != scalar_n_dofs; k++)
465+
for (unsigned int n = 0; n != lambda_n_dofs; n++)
466+
{
467+
Ke(k + vector_n_dofs, n + vector_n_dofs + scalar_n_dofs) += JxW[qp]*scalar_phi[k][qp];
468+
}
469+
470+
// A double loop to integrate the Lagrange dof (m) against the
471+
// scalar trial functions (l).
472+
for (unsigned int m = 0; m != lambda_n_dofs; m++)
459473
for (unsigned int l = 0; l != scalar_n_dofs; l++)
460474
{
461-
Ke(k + vector_n_dofs, l + vector_n_dofs) += JxW[qp]*scalar_phi[k][qp]*scalar_phi[l][qp];
475+
Ke(m + vector_n_dofs + scalar_n_dofs, l + vector_n_dofs) += JxW[qp]*scalar_phi[l][qp];
462476
}
463477

464-
// Loop to integrate the scalar test functions (k) against the
465-
// exact solution for the scalar variable.
466-
for (unsigned int k = 0; k != scalar_n_dofs; k++)
478+
// Loop to integrate the exact solution for the scalar variable.
479+
for (unsigned int m = 0; m != lambda_n_dofs; m++)
467480
{
468-
Fe(k + vector_n_dofs) += JxW[qp]*scalar_phi[k][qp]*scalar_value;
481+
Fe(m + vector_n_dofs + scalar_n_dofs) += JxW[qp]*scalar_value;
469482
}
470483
}
471484
}

0 commit comments

Comments
 (0)