GetFEM  5.5
getfem_HHO.cc
1 /*===========================================================================
2 
3  Copyright (C) 2019-2026 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program. If not, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
21 
22 #include "getfem/getfem_HHO.h"
23 
24 
25 namespace getfem {
26 
27  THREAD_SAFE_STATIC bgeot::geotrans_precomp_pool HHO_pgp_pool;
28  THREAD_SAFE_STATIC fem_precomp_pool HHO_pfp_pool;
29 
30  // To be optimized:
31  // - The fact that (when pf2->target_dim() = 1) the
32  // problem can be solved componentwise can be more exploited in
33  // avoiding the computation of the whole matrix M2.
34  // - The vectorization can be avoided in most cases
35  // - Avoid some multiple ctx
36 
37 
38  class _HHO_reconstructed_gradient
39  : public virtual_elementary_transformation {
40 
41  public:
42 
43  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
44  size_type cv, base_matrix &M) const {
45 
46  // The reconstructed Gradient "G" is described on mf2 and computed by
47  // the formula on the element T :
48  // \int_T G.w = \int_T Grad(v_T).w + \int_{dT}(v_{dT} - v_T).(w.n)
49  // where "w" is the test function arbitrary in mf2, "v_T" is the field
50  // inside the element whose gradient is to be reconstructed,
51  // "v_{dT}" is the field on the boundary of T and "n" is the outward
52  // unit normal.
53 
54  // Obtaining the fem descriptors
55  pfem pf1 = mf1.fem_of_element(cv);
56  pfem pf2 = mf2.fem_of_element(cv);
57  pfem pfi = interior_fem_of_hho_method(pf1);
58 
59  size_type degree = std::max(pf1->estimated_degree(),
60  pf2->estimated_degree());
61  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
62  papprox_integration pim
63  = classical_approx_im(pgt, dim_type(2*degree))->approx_method();
64 
65  base_matrix G;
66  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
67 
68  bgeot::pgeotrans_precomp pgp
69  = HHO_pgp_pool(pgt, pim->pintegration_points());
70  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
71  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
72  pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
73 
74  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
75  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
76  fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
77 
78  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
79  base_vector un(N);
80  size_type qmult1 = Q / pf1->target_dim();
81  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
82  size_type qmult2 = Q*N / pf2->target_dim();
83  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
84  size_type qmulti = Q / pfi->target_dim();
85  size_type ndofi = pfi->nb_dof(cv) * qmulti;
86 
87  base_tensor t1, t2, ti, tv1;
88  base_matrix tv2, tv1p, tvi;
89  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
90  base_matrix aux2(ndof2, ndof2);
91 
92  // Integrals on the element : \int_T G.w (M2) and \int_T Grad(v_T).w (M1)
93  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
94  ctx1.set_ii(ipt); ctx2.set_ii(ipt);
95  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
96 
97  ctx1.grad_base_value(t1);
98  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
99 
100  ctx2.base_value(t2);
101  vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q*N);
102 
103  gmm::mult(tv2, gmm::transposed(tv2), aux2);
104  gmm::add(gmm::scaled(aux2, coeff), M2);
105 
106  for (size_type i = 0; i < ndof1; ++i) // To be optimized
107  for (size_type j = 0; j < ndof2; ++j)
108  for (size_type k = 0; k < Q*N; ++k)
109  M1(j, i) += coeff * tv1.as_vector()[i+k*ndof1] * tv2(j, k);
110  }
111 
112  // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(w.n) (M1)
113  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
114  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
115  size_type first_ind = pim->ind_first_point_on_face(ifc);
116  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
117  ctx1.set_ii(first_ind+ipt);
118  ctx2.set_ii(first_ind+ipt);
119  ctxi.set_ii(first_ind+ipt);
120  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
121  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
122 
123  ctx2.base_value(t2);
124  vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q*N);
125 
126  ctx1.base_value(t1);
127  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
128 
129  ctxi.base_value(ti);
130  vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), Q);
131 
132  for (size_type i = 0; i < ndof1; ++i) // To be optimized
133  for (size_type j = 0; j < ndof2; ++j)
134  for (size_type k1 = 0; k1 < Q; ++k1) {
135  scalar_type b(0), a = coeff *
136  (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
137  for (size_type k2 = 0; k2 < N; ++k2)
138  b += a * tv2(j, k1 + k2*Q) * un[k2];
139  M1(j, i) += b;
140  }
141  }
142  }
143 
144  if (pf2->target_dim() == 1) {
145  gmm::sub_slice I(0, ndof2/(N*Q), N*Q);
146  gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
147  for (size_type i = 0; i < N*Q; ++i) {
148  gmm::sub_slice I2(i, ndof2/(N*Q), N*Q);
149  gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
150  }
151  } else { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
152 
153  gmm::mult(M2inv, M1, M);
154  gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
155  }
156  };
157 
158  void add_HHO_reconstructed_gradient(model &md, std::string name) {
159  pelementary_transformation
160  p = std::make_shared<_HHO_reconstructed_gradient>();
161  md.add_elementary_transformation(name, p);
162  }
163 
164 
165  class _HHO_reconstructed_sym_gradient
166  : public virtual_elementary_transformation {
167 
168  public:
169 
170  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
171  size_type cv, base_matrix &M) const {
172 
173  // The reconstructed symmetric Gradient "G" is described on mf2 and
174  // computed by the formula on the element T :
175  // \int_T G:w = (1/2)*\int_T 0.5*Grad(v_T):(w+w^T)
176  // + (1/2)*\int_{dT}(v_{dT} - v_T).((w+w^T).n)
177  // where "w" is the test function arbitrary in mf2, "v_T" is the field
178  // inside the element whose gradient is to be reconstructed,
179  // "v_{dT}" is the field on the boundary of T and "n" is the outward
180  // unit normal.
181 
182  // Obtaining the fem descriptors
183  pfem pf1 = mf1.fem_of_element(cv);
184  pfem pf2 = mf2.fem_of_element(cv);
185  pfem pfi = interior_fem_of_hho_method(pf1);
186 
187  size_type degree = std::max(pf1->estimated_degree(),
188  pf2->estimated_degree());
189  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
190  papprox_integration pim
191  = classical_approx_im(pgt, dim_type(2*degree))->approx_method();
192 
193  base_matrix G;
194  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
195 
196  bgeot::pgeotrans_precomp pgp
197  = HHO_pgp_pool(pgt, pim->pintegration_points());
198  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
199  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
200  pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
201 
202  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
203  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
204  fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
205 
206  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
207  GMM_ASSERT1(Q == N, "This transformation works only for vector fields "
208  "having the same dimension as the domain");
209  base_vector un(N);
210  size_type qmult1 = N / pf1->target_dim();
211  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
212  size_type qmult2 = N*N / pf2->target_dim();
213  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
214  size_type qmulti = N / pfi->target_dim();
215  size_type ndofi = pfi->nb_dof(cv) * qmulti;
216 
217 
218  base_tensor t1, t2, ti, tv1;
219  base_matrix tv2, tv1p, tvi;
220  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
221  base_matrix aux2(ndof2, ndof2);
222 
223  // Integrals on the element : \int_T G:w (M2)
224  // and (1/2)*\int_T 0.5*Grad(v_T):(w+w^T)
225  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
226  ctx1.set_ii(ipt); ctx2.set_ii(ipt);
227  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
228 
229  ctx1.grad_base_value(t1);
230  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
231 
232  ctx2.base_value(t2);
233  vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N*N);
234 
235  gmm::mult(tv2, gmm::transposed(tv2), aux2);
236  gmm::add(gmm::scaled(aux2, coeff), M2);
237 
238  for (size_type i = 0; i < ndof1; ++i) // To be optimized
239  for (size_type j = 0; j < ndof2; ++j)
240  for (size_type k1 = 0; k1 < N; ++k1)
241  for (size_type k2 = 0; k2 < N; ++k2)
242  M1(j, i) += coeff * tv1.as_vector()[i+(k1+k2*N)*ndof1]
243  * 0.5 * (tv2(j, k1+k2*N) + tv2(j, k2+k1*N));
244  }
245 
246  // Integrals on the faces : (1/2)*\int_{dT}(v_{dT} - v_T).((w+w^T).n) (M1)
247  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
248  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
249  size_type first_ind = pim->ind_first_point_on_face(ifc);
250  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
251  ctx1.set_ii(first_ind+ipt);
252  ctx2.set_ii(first_ind+ipt);
253  ctxi.set_ii(first_ind+ipt);
254  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
255  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
256 
257  ctx2.base_value(t2);
258  vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N*N);
259 
260  ctx1.base_value(t1);
261  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), N);
262 
263  ctxi.base_value(ti);
264  vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), N);
265 
266  for (size_type i = 0; i < ndof1; ++i) // To be optimized
267  for (size_type j = 0; j < ndof2; ++j)
268  for (size_type k1 = 0; k1 < N; ++k1) {
269  scalar_type b(0), a = coeff *
270  (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
271  for (size_type k2 = 0; k2 < N; ++k2)
272  b += a*0.5*(tv2(j, k1 + k2*N) + tv2(j, k2 + k1*N)) * un[k2];
273  M1(j, i) += b;
274  }
275  }
276 
277  }
278 
279  if (pf2->target_dim() == 1) {
280  gmm::sub_slice I(0, ndof2/(N*Q), N*Q);
281  gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
282  for (size_type i = 0; i < N*Q; ++i) {
283  gmm::sub_slice I2(i, ndof2/(N*Q), N*Q);
284  gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
285  }
286  } else
287  { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
288 
289  gmm::mult(M2inv, M1, M);
290  gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
291  }
292  };
293 
294  void add_HHO_reconstructed_symmetrized_gradient(model &md, std::string name) {
295  pelementary_transformation
296  p = std::make_shared<_HHO_reconstructed_sym_gradient>();
297  md.add_elementary_transformation(name, p);
298  }
299 
300 
301 
302  class _HHO_reconstructed_value
303  : public virtual_elementary_transformation {
304 
305  public:
306 
307  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
308  size_type cv, base_matrix &M) const {
309  // The reconstructed variable "D" is described on mf2 and computed by
310  // the formula on the element T :
311  // \int_T Grad(D).Grad(w) = \int_T Grad(v_T).Grad(w)
312  // + \int_{dT}(v_{dT} - v_T).(Grad(w).n)
313  // with the constraint
314  // \int_T D = \int_T v_T
315  // where "w" is the test function arbitrary in mf2, "v_T" is the field
316  // inside the element whose gradient is to be reconstructed,
317  // "v_{dT}" is the field on the boundary of T and "n" is the outward
318  // unit normal.
319 
320  // Obtaining the fem descriptors
321  pfem pf1 = mf1.fem_of_element(cv);
322  pfem pf2 = mf2.fem_of_element(cv);
323  pfem pfi = interior_fem_of_hho_method(pf1);
324 
325  size_type degree = std::max(pf1->estimated_degree(),
326  pf2->estimated_degree());
327  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
328  papprox_integration pim
329  = classical_approx_im(pgt, dim_type(2*degree))->approx_method();
330 
331  base_matrix G;
332  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
333 
334  bgeot::pgeotrans_precomp pgp
335  = HHO_pgp_pool(pgt, pim->pintegration_points());
336  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
337  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
338  pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
339 
340  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
341  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
342  fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
343 
344  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
345  base_vector un(N);
346  size_type qmult1 = Q / pf1->target_dim();
347  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
348  size_type qmult2 = Q / pf2->target_dim();
349  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
350  size_type qmulti = Q / pfi->target_dim();
351  size_type ndofi = pfi->nb_dof(cv) * qmulti;
352 
353 
354  base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
355  base_matrix tv1p, tv2p, tvi;
356  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
357  base_matrix M3(Q, ndof1), M4(Q, ndof2);
358  base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
359  scalar_type area(0);
360 
361  // Integrals on the element : \int_T Grad(D).Grad(w) (M2)
362  // \int_T Grad(v_T).Grad(w) (M1)
363  // \int_T D (M4) and \int_T v_T (M3)
364  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
365  ctx1.set_ii(ipt); ctx2.set_ii(ipt);
366  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
367  area += coeff;
368 
369  ctx1.grad_base_value(t1);
370  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
371 
372  ctx1.base_value(t1p);
373  vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), Q);
374 
375  ctx2.grad_base_value(t2);
376  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
377 
378  ctx2.base_value(t2p);
379  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
380 
381  for (size_type i = 0; i < ndof2; ++i) // To be optimized
382  for (size_type j = 0; j < ndof2; ++j)
383  for (size_type k = 0; k < Q*N; ++k)
384  M2(j, i) += coeff * tv2.as_vector()[i+k*ndof2]
385  * tv2.as_vector()[j+k*ndof2];
386 
387  for (size_type i = 0; i < ndof2; ++i)
388  for (size_type k = 0; k < Q; ++k)
389  M4(k, i) += coeff * tv2p(i, k);
390 
391  for (size_type i = 0; i < ndof1; ++i) // To be optimized
392  for (size_type j = 0; j < ndof2; ++j)
393  for (size_type k = 0; k < Q*N; ++k)
394  M1(j, i) += coeff * tv1.as_vector()[i+k*ndof1]
395  * tv2.as_vector()[j+k*ndof2];
396 
397  for (size_type i = 0; i < ndof1; ++i)
398  for (size_type k = 0; k < Q; ++k)
399  M3(k, i) += coeff * tv1p(i, k);
400 
401  }
402 
403  // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
404  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
405  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
406  size_type first_ind = pim->ind_first_point_on_face(ifc);
407  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
408  ctx1.set_ii(first_ind+ipt);
409  ctx2.set_ii(first_ind+ipt);
410  ctxi.set_ii(first_ind+ipt);
411  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
412  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
413 
414  ctx2.grad_base_value(t2);
415  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
416 
417  ctx1.base_value(t1);
418  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
419 
420  ctxi.base_value(ti);
421  vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), Q);
422 
423  for (size_type i = 0; i < ndof1; ++i) // To be optimized
424  for (size_type j = 0; j < ndof2; ++j)
425  for (size_type k1 = 0; k1 < Q; ++k1) {
426  scalar_type b(0), a = coeff *
427  (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
428  for (size_type k2 = 0; k2 < N; ++k2)
429  b += a * tv2.as_vector()[j+(k1+k2*Q)*ndof2] * un[k2];
430  M1(j, i) += b;
431  }
432  }
433 
434  }
435 
436  // Add the constraint with penalization
437  scalar_type coeff_p = pow(area, -1. - 2./scalar_type(N));
438  gmm::mult(gmm::transposed(M4), M4, aux2);
439  gmm::add (gmm::scaled(aux2, coeff_p), M2);
440  gmm::mult(gmm::transposed(M4), M3, aux1);
441  gmm::add (gmm::scaled(aux1, coeff_p), M1);
442 
443  if (pf2->target_dim() == 1 && Q > 1) {
444  gmm::sub_slice I(0, ndof2/Q, Q);
445  gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
446  for (size_type i = 0; i < Q; ++i) {
447  gmm::sub_slice I2(i, ndof2/Q, Q);
448  gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
449  }
450  } else
451  { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
452 
453  gmm::mult(M2inv, M1, M);
454  gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
455  }
456  };
457 
458  void add_HHO_reconstructed_value(model &md, std::string name) {
459  pelementary_transformation
460  p = std::make_shared<_HHO_reconstructed_value>();
461  md.add_elementary_transformation(name, p);
462  }
463 
464 
465  class _HHO_reconstructed_sym_value
466  : public virtual_elementary_transformation {
467 
468  public:
469 
470  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
471  size_type cv, base_matrix &M) const {
472  // The reconstructed variable "D" is described on mf2 and computed by
473  // the formula on the element T :
474  // \int_T Sym(Grad(D)).Grad(w) = \int_T Sym(Grad(v_T)).Grad(w)
475  // + \int_{dT}(v_{dT} - v_T).(Sym(Grad(w)).n)
476  // with the constraints
477  // \int_T D = \int_T v_T
478  // \int_T Skew(Grad(D)) = 0.5\int_{dT}(n x v_{dT} - v_{dT} x n)
479  // where "w" is the test function arbitrary in mf2, "v_T" is the field
480  // inside the element whose gradient is to be reconstructed,
481  // "v_{dT}" is the field on the boundary of T and "n" is the outward
482  // unit normal.
483 
484  // Obtaining the fem descriptors
485  pfem pf1 = mf1.fem_of_element(cv);
486  pfem pf2 = mf2.fem_of_element(cv);
487  pfem pfi = interior_fem_of_hho_method(pf1);
488 
489  size_type degree = std::max(pf1->estimated_degree(),
490  pf2->estimated_degree());
491  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
492  papprox_integration pim
493  = classical_approx_im(pgt, dim_type(2*degree))->approx_method();
494 
495  base_matrix G;
496  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
497 
498  bgeot::pgeotrans_precomp pgp
499  = HHO_pgp_pool(pgt, pim->pintegration_points());
500  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
501  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
502  pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
503 
504  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
505  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
506  fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
507 
508  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
509  GMM_ASSERT1(Q == N, "This transformation works only for vector fields "
510  "having the same dimension as the domain");
511  base_vector un(N);
512  size_type qmult1 = N / pf1->target_dim();
513  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
514  size_type qmult2 = N / pf2->target_dim();
515  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
516  size_type qmulti = N / pfi->target_dim();
517  size_type ndofi = pfi->nb_dof(cv) * qmulti;
518 
519 
520  base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
521  base_matrix tv1p, tv2p, tvi;
522  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);;
523  base_matrix M3(N, ndof1), M4(N, ndof2);
524  base_matrix M5(N*N, ndof1), M6(N*N, ndof2);
525  base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
526  scalar_type area(0);
527 
528  // Integrals on the element : \int_T Sym(Grad(D)).Grad(w) (M2)
529  // \int_T Sym(Grad(v_T)).Grad(w) (M1)
530  // \int_T D (M4) and \int_T v_T (M3)
531  // \int_T Skew(Grad(D)) (M6)
532  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
533  ctx1.set_ii(ipt); ctx2.set_ii(ipt);
534  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
535  area += coeff;
536 
537  ctx1.grad_base_value(t1);
538  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
539 
540  ctx1.base_value(t1p);
541  vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), N);
542 
543  ctx2.grad_base_value(t2);
544  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
545 
546  ctx2.base_value(t2p);
547  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), N);
548 
549  for (size_type i = 0; i < ndof2; ++i) // To be optimized
550  for (size_type j = 0; j < ndof2; ++j)
551  for (size_type k1 = 0; k1 < N; ++k1)
552  for (size_type k2 = 0; k2 < N; ++k2)
553  M2(j, i) += coeff * tv2.as_vector()[i+(k1+k2*N)*ndof2]
554  * 0.5 * (tv2.as_vector()[j+(k1+k2*N)*ndof2] +
555  tv2.as_vector()[j+(k2+k1*N)*ndof2]);
556 
557  for (size_type i = 0; i < ndof2; ++i)
558  for (size_type k = 0; k < N; ++k)
559  M4(k, i) += coeff * tv2p(i, k);
560 
561  for (size_type i = 0; i < ndof2; ++i)
562  for (size_type k1 = 0; k1 < N; ++k1)
563  for (size_type k2 = 0; k2 < N; ++k2)
564  M6(k1+k2*N, i) += 0.5*coeff*(tv2.as_vector()[i+(k1+k2*N)*ndof2] -
565  tv2.as_vector()[i+(k2+k1*N)*ndof2]);
566 
567  for (size_type i = 0; i < ndof1; ++i) // To be optimized
568  for (size_type j = 0; j < ndof2; ++j)
569  for (size_type k1 = 0; k1 < N; ++k1)
570  for (size_type k2 = 0; k2 < N; ++k2)
571  M1(j, i) += coeff * tv1.as_vector()[i+(k1+k2*N)*ndof1]
572  * 0.5 * (tv2.as_vector()[j+(k1+k2*N)*ndof2] +
573  tv2.as_vector()[j+(k2+k1*N)*ndof2]);
574 
575  for (size_type i = 0; i < ndof1; ++i)
576  for (size_type k = 0; k < N; ++k)
577  M3(k, i) += coeff * tv1p(i, k);
578 
579  }
580 
581  // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Sym(Grad(w)).n) (M1)
582  // \int_{dT} n x v_{dT} - v_{dT} x n (M5)
583  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
584  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
585  size_type first_ind = pim->ind_first_point_on_face(ifc);
586  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
587  ctx1.set_ii(first_ind+ipt);
588  ctx2.set_ii(first_ind+ipt);
589  ctxi.set_ii(first_ind+ipt);
590  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
591  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
592 
593  ctx2.grad_base_value(t2);
594  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
595 
596  ctx1.base_value(t1);
597  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), N);
598 
599  ctxi.base_value(ti);
600  vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), N);
601 
602  for (size_type i = 0; i < ndof1; ++i) // To be optimized
603  for (size_type j = 0; j < ndof2; ++j)
604  for (size_type k1 = 0; k1 < N; ++k1) {
605  scalar_type b(0), a = coeff *
606  (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
607  for (size_type k2 = 0; k2 < N; ++k2)
608  b += a * 0.5 * (tv2.as_vector()[j+(k1+k2*N)*ndof2] +
609  tv2.as_vector()[j+(k2+k1*N)*ndof2]) * un[k2];
610  M1(j, i) += b;
611  }
612 
613  for (size_type i = 0; i < ndof1; ++i)
614  for (size_type k1 = 0; k1 < N; ++k1)
615  for (size_type k2 = 0; k2 < N; ++k2)
616  M5(k1+k2*N, i) += 0.5 * coeff * (tv1p(i, k1) * un[k2] -
617  tv1p(i, k2) * un[k1]);
618  }
619  }
620 
621  // Add the constraint with penalization
622  scalar_type coeff_p1 = pow(area, -1. - 2./scalar_type(N));
623  scalar_type coeff_p2 = pow(area, -1. - 1./scalar_type(N));
624  gmm::mult(gmm::transposed(M4), M4, aux2);
625  gmm::add (gmm::scaled(aux2, coeff_p1), M2);
626  gmm::mult(gmm::transposed(M6), M6, aux2);
627  gmm::add (gmm::scaled(aux2, coeff_p2), M2);
628  gmm::mult(gmm::transposed(M4), M3, aux1);
629  gmm::add (gmm::scaled(aux1, coeff_p1), M1);
630  gmm::mult(gmm::transposed(M6), M5, aux1);
631  gmm::add (gmm::scaled(aux1, coeff_p2), M1);
632 
633  gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv);
634 
635  gmm::mult(M2inv, M1, M);
636  gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
637  }
638  };
639 
640  void add_HHO_reconstructed_symmetrized_value(model &md, std::string name) {
641  pelementary_transformation
642  p = std::make_shared<_HHO_reconstructed_sym_value>();
643  md.add_elementary_transformation(name, p);
644  }
645 
646 #if 0 // Old single mef version
647 
648 class _HHO_stabilization
649  : public virtual_elementary_transformation {
650 
651  public:
652 
653  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
654  size_type cv, base_matrix &M) const {
655  // The reconstructed variable "S" is described on mf2 and computed by
656  // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)))
657  // where P__{\dT} et P_T are L2 projections on the boundary and on the
658  // interior of T on the corresponding discrete spaces.
659  // Note that P_{\dT}(v_{dT}) = v_{dT} and P_T(v_T) = v_T and D is
660  // the reconstructed value on P^{k+1} given by the formula:
661  // \int_T Grad(D).Grad(w) = \int_T Grad(v_T).Grad(w)
662  // + \int_{dT}(v_{dT} - v_T).(Grad(w).n)
663  // with the constraint
664  // \int_T D = \int_T v_T
665  // where "w" is the test function arbitrary in mf2, "v_T" is the field
666  // inside the element whose gradient is to be reconstructed,
667  // "v_{dT}" is the field on the boundary of T and "n" is the outward
668  // unit normal.
669  // The implemented formula is
670  // S(v) = v_{dT} - P_{\dT}D(v) - P_{\dT}(v_T) + P_{\dT}(P_T(D(v)))
671  // by the mean of the projection matrix from P^{k+1} to the original space
672  // and the projection matrix from interior space to the boundary space.
673  // As it is built, S(v) is zero on interior dofs.
674 
675  GMM_ASSERT1(&mf1 == &mf2, "The HHO stabilization transformation is "
676  "only defined on the HHO space to itself");
677 
678  // Obtaining the fem descriptors
679  pfem pf1 = mf1.fem_of_element(cv);
680  short_type degree = pf1->estimated_degree();
681  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
682  pfem pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be
683  // changed for an interior PK method
684  pfem pfi = interior_fem_of_hho_method(pf1);
685 
686  papprox_integration pim
687  = classical_approx_im(pgt, dim_type(2*degree+2))->approx_method();
688 
689  base_matrix G;
690  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
691 
692  bgeot::pgeotrans_precomp pgp
693  = HHO_pgp_pool(pgt, pim->pintegration_points());
694  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
695  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
696  pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
697 
698  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
699  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
700  fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
701 
702  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
703  base_vector un(N);
704  size_type qmult1 = Q / pf1->target_dim();
705  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
706  size_type qmult2 = Q / pf2->target_dim();
707  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
708  size_type qmulti = Q / pfi->target_dim();
709  size_type ndofi = pfi->nb_dof(cv) * qmulti;
710 
711 
712  base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
713  base_matrix tv1p, tv2p, tvi;
714  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
715  base_matrix M3(Q, ndof1), M4(Q, ndof2);
716  base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
717  base_matrix M7(ndof1, ndof1), M7inv(ndof1, ndof1), M8(ndof1, ndof2);
718  base_matrix M9(ndof1, ndof1), MD(ndof2, ndof1);
719  scalar_type area(0);
720 
721  // Integrals on the element : \int_T Grad(D).Grad(w) (M2)
722  // \int_T Grad(v_T).Grad(w) (M1)
723  // \int_T D (M4) and \int_T v_T (M3)
724  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
725  ctx1.set_ii(ipt); ctx2.set_ii(ipt);
726  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
727  area += coeff;
728 
729  ctx1.grad_base_value(t1);
730  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
731 
732  ctx1.base_value(t1p);
733  vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), Q);
734 
735  ctx2.grad_base_value(t2);
736  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
737 
738  ctx2.base_value(t2p);
739  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
740 
741  for (size_type i = 0; i < ndof2; ++i) // To be optimized
742  for (size_type j = 0; j < ndof2; ++j)
743  for (size_type k = 0; k < Q*N; ++k)
744  M2(j, i) += coeff * tv2.as_vector()[i+k*ndof2]
745  * tv2.as_vector()[j+k*ndof2];
746 
747  for (size_type i = 0; i < ndof2; ++i)
748  for (size_type k = 0; k < Q; ++k)
749  M4(k, i) += coeff * tv2p(i, k);
750 
751  for (size_type i = 0; i < ndof1; ++i) // To be optimized
752  for (size_type j = 0; j < ndof2; ++j)
753  for (size_type k = 0; k < Q*N; ++k)
754  M1(j, i) += coeff * tv1.as_vector()[i+k*ndof1]
755  * tv2.as_vector()[j+k*ndof2];
756 
757  for (size_type i = 0; i < ndof1; ++i)
758  for (size_type k = 0; k < Q; ++k)
759  M3(k, i) += coeff * tv1p(i, k);
760 
761  for (size_type i = 0; i < ndof1; ++i) // To be optimized
762  for (size_type j = 0; j < ndof1; ++j)
763  for (size_type k = 0; k < Q; ++k)
764  M7(i, j) += coeff * tv1p(i, k) * tv1p(j, k);
765 
766  for (size_type i = 0; i < ndof1; ++i) // To be optimized
767  for (size_type j = 0; j < ndof2; ++j)
768  for (size_type k = 0; k < Q; ++k)
769  M8(i, j) += coeff * tv1p(i, k) * tv2p(j, k);
770 
771  }
772 
773  // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
774  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
775  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
776  size_type first_ind = pim->ind_first_point_on_face(ifc);
777  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
778  ctx1.set_ii(first_ind+ipt);
779  ctx2.set_ii(first_ind+ipt);
780  ctxi.set_ii(first_ind+ipt);
781  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
782  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
783  scalar_type normun = gmm::vect_norm2(un);
784 
785  ctx2.grad_base_value(t2);
786  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
787 
788  ctx2.base_value(t2p);
789  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
790 
791  ctx1.base_value(t1);
792  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
793 
794  ctxi.base_value(ti);
795  vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), Q);
796 
797 
798  for (size_type i = 0; i < ndof1; ++i) // To be optimized
799  for (size_type j = 0; j < ndof2; ++j)
800  for (size_type k1 = 0; k1 < Q; ++k1) {
801  scalar_type b(0), a = coeff *
802  (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
803  for (size_type k2 = 0; k2 < N; ++k2)
804  b += a * tv2.as_vector()[j+(k1 + k2*Q)*ndof2] * un[k2];
805  M1(j, i) += b;
806  }
807 
808  for (size_type i = 0; i < ndof1; ++i) // To be optimized
809  for (size_type j = 0; j < ndof1; ++j)
810  for (size_type k = 0; k < Q; ++k)
811  M7(i, j) += coeff * normun * tv1p(i,k) * tv1p(j, k);
812 
813  for (size_type i = 0; i < ndof1; ++i) // To be optimized
814  for (size_type j = 0; j < ndof2; ++j)
815  for (size_type k = 0; k < Q; ++k)
816  M8(i, j) += coeff * normun * tv1p(i,k) * tv2p(j, k);
817 
818  for (size_type i = 0; i < ndof1; ++i) // To be optimized
819  for (size_type j = 0; j < ndofi; ++j)
820  for (size_type k = 0; k < Q; ++k)
821  M9(i, j) += coeff * normun * tv1p(i,k) * tvi(j, k);
822  }
823  }
824 
825  // Add the constraint with penalization
826  scalar_type coeff_p = pow(area, -1. - 2./scalar_type(N));
827  gmm::mult(gmm::transposed(M4), M4, aux2);
828  gmm::add (gmm::scaled(aux2, coeff_p), M2);
829  gmm::mult(gmm::transposed(M4), M3, aux1);
830  gmm::add (gmm::scaled(aux1, coeff_p), M1);
831 
832  if (pf2->target_dim() == 1 && Q > 1) {
833  gmm::sub_slice I(0, ndof2/Q, Q);
834  gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
835  for (size_type i = 0; i < Q; ++i) {
836  gmm::sub_slice I2(i, ndof2/Q, Q);
837  gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
838  }
839  } else
840  { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
841 
842  if (pf1->target_dim() == 1 && Q > 1) {
843  gmm::sub_slice I(0, ndof1/Q, Q);
844  gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
845  for (size_type i = 0; i < Q; ++i) {
846  gmm::sub_slice I2(i, ndof1/Q, Q);
847  gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
848  }
849  } else
850  { gmm::copy(M7, M7inv); gmm::lu_inverse(M7inv); }
851 
852  gmm::mult(M2inv, M1, MD);
853  gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
854 
855  // S = (I - inv(M7)*M9)(I - inv(M7)*M8*MD)
856  base_matrix MPB(ndof1, ndof1);
857  gmm::mult(M7inv, M9, MPB);
858  gmm::copy(gmm::identity_matrix(), M9);
859  gmm::add(gmm::scaled(MPB, scalar_type(-1)), M9);
860 
861  base_matrix MPC(ndof1, ndof1), MPD(ndof1, ndof1);
862  gmm::mult(M8, MD, MPC);
863  gmm::mult(M7inv, MPC, MPD);
864  gmm::copy(gmm::identity_matrix(), M7);
865  gmm::add(gmm::scaled(MPD, scalar_type(-1)), M7);
866 
867  gmm::mult(M9, M7, M);
868  gmm::clean(M, 1E-13);
869  }
870  };
871 
872  void add_HHO_stabilization(model &md, std::string name) {
873  pelementary_transformation
874  p = std::make_shared<_HHO_stabilization>();
875  md.add_elementary_transformation(name, p);
876  }
877 
878 
879 #else
880 
881 
882  class _HHO_stabilization
883  : public virtual_elementary_transformation {
884 
885  public:
886 
887  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
888  size_type cv, base_matrix &M) const {
889  // The reconstructed variable "S" is described on mf2 and computed by
890  // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)))
891  // where P_{\dT} et P_T are L2 projections on the boundary and on the
892  // interior of T on the corresponding discrete spaces.
893  // Note that P_{\dT}(v_{dT}) = v_{dT} and P_T(v_T) = v_T and D is
894  // the reconstructed value on P^{k+1} given by the formula:
895  // \int_T Grad(D).Grad(w) = \int_T Grad(v_T).Grad(w)
896  // + \int_{dT}(v_{dT} - v_T).(Grad(w).n)
897  // with the constraint
898  // \int_T D = \int_T v_T
899  // where "w" is the test function arbitrary in mf2, "v_T" is the field
900  // inside the element whose gradient is to be reconstructed,
901  // "v_{dT}" is the field on the boundary of T and "n" is the outward
902  // unit normal.
903  // The implemented formula is
904  // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)) )
905  // by the mean of the projection matrix from P^{k+1} to the target space
906  // and the projection matrix from interior space to the boundary space.
907  // As it is built, S(v) is zero on interior dofs.
908 
909  // Obtaining the fem descriptors
910  pfem pf1 = mf1.fem_of_element(cv);
911  short_type degree = pf1->estimated_degree();
912  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
913  pfem pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be
914  // changed for an interior PK method
915  pfem pf3 = mf2.fem_of_element(cv);
916  pfem pf1i = interior_fem_of_hho_method(pf1);
917  pfem pf3i = interior_fem_of_hho_method(pf3);
918 
919  papprox_integration pim
920  = classical_approx_im(pgt, dim_type(2*degree+2))->approx_method();
921 
922  base_matrix G;
923  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
924 
925  bgeot::pgeotrans_precomp pgp
926  = HHO_pgp_pool(pgt, pim->pintegration_points());
927  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
928  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
929  pfem_precomp pfp3 = HHO_pfp_pool(pf3, pim->pintegration_points());
930  pfem_precomp pfp1i = HHO_pfp_pool(pf1i, pim->pintegration_points());
931  pfem_precomp pfp3i = HHO_pfp_pool(pf3i, pim->pintegration_points());
932 
933  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
934  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
935  fem_interpolation_context ctx3(pgp, pfp3, 0, G, cv);
936  fem_interpolation_context ctx1i(pgp, pfp1i, 0, G, cv);
937  fem_interpolation_context ctx3i(pgp, pfp3i, 0, G, cv);
938 
939  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
940  base_vector un(N);
941  size_type qmult1 = Q / pf1->target_dim();
942  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
943  size_type qmult2 = Q / pf2->target_dim();
944  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
945  size_type qmult3 = Q / pf3->target_dim();
946  size_type ndof3 = pf3->nb_dof(cv) * qmult3;
947  size_type qmult1i = Q / pf1i->target_dim();
948  size_type ndof1i = pf1i->nb_dof(cv) * qmult1i;
949  size_type qmult3i = Q / pf3i->target_dim();
950  size_type ndof3i = pf3i->nb_dof(cv) * qmult3i;
951 
952 
953  base_tensor t1, t2, t3, t1i, t3i, tv1, tv2, t1p, t2p;
954  base_matrix tv1p, tv2p, tv3p, tv1i, tv3i;
955  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
956  base_matrix M3(Q, ndof1), M4(Q, ndof2);
957  base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
958  base_matrix M7(ndof3, ndof3), M7inv(ndof3, ndof3), M8(ndof3, ndof2);
959  base_matrix M9(ndof3, ndof1), M10(ndof3, ndof3), MD(ndof2, ndof1);
960  scalar_type area(0);
961 
962  // Integrals on the element : \int_T Grad(D).Grad(w) (M2)
963  // \int_T Grad(v_T).Grad(w) (M1)
964  // \int_T D (M4) and \int_T v_T (M3)
965  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
966  ctx1.set_ii(ipt); ctx2.set_ii(ipt); ctx3.set_ii(ipt);
967  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
968  area += coeff;
969 
970  ctx1.grad_base_value(t1);
971  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
972 
973  ctx1.base_value(t1p);
974  vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), Q);
975 
976  ctx2.grad_base_value(t2);
977  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
978 
979  ctx2.base_value(t2p);
980  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
981 
982  ctx3.base_value(t3);
983  vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
984 
985  for (size_type i = 0; i < ndof2; ++i) // To be optimized
986  for (size_type j = 0; j < ndof2; ++j)
987  for (size_type k = 0; k < Q*N; ++k)
988  M2(j, i) += coeff * tv2.as_vector()[i+k*ndof2]
989  * tv2.as_vector()[j+k*ndof2];
990 
991  for (size_type i = 0; i < ndof2; ++i)
992  for (size_type k = 0; k < Q; ++k)
993  M4(k, i) += coeff * tv2p(i, k);
994 
995  for (size_type i = 0; i < ndof1; ++i) // To be optimized
996  for (size_type j = 0; j < ndof2; ++j)
997  for (size_type k = 0; k < Q*N; ++k)
998  M1(j, i) += coeff * tv1.as_vector()[i+k*ndof1]
999  * tv2.as_vector()[j+k*ndof2];
1000 
1001  for (size_type i = 0; i < ndof1; ++i)
1002  for (size_type k = 0; k < Q; ++k)
1003  M3(k, i) += coeff * tv1p(i, k);
1004 
1005  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1006  for (size_type j = 0; j < ndof3; ++j)
1007  for (size_type k = 0; k < Q; ++k)
1008  M7(i, j) += coeff * tv3p(i, k) * tv3p(j, k);
1009 
1010  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1011  for (size_type j = 0; j < ndof2; ++j)
1012  for (size_type k = 0; k < Q; ++k)
1013  M8(i, j) += coeff * tv3p(i, k) * tv2p(j, k);
1014 
1015  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1016  for (size_type j = 0; j < ndof1; ++j)
1017  for (size_type k = 0; k < Q; ++k)
1018  M9(i, j) += coeff * tv3p(i, k) * tv1p(j, k);
1019  }
1020 
1021  // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
1022  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
1023  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctx3.set_face_num(ifc);
1024  ctx1i.set_face_num(ifc); ctx3i.set_face_num(ifc);
1025  size_type first_ind = pim->ind_first_point_on_face(ifc);
1026  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
1027  ctx1.set_ii(first_ind+ipt); ctx2.set_ii(first_ind+ipt);
1028  ctx3.set_ii(first_ind+ipt); ctx1i.set_ii(first_ind+ipt);
1029  ctx3i.set_ii(first_ind+ipt);
1030  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
1031  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
1032  scalar_type normun = gmm::vect_norm2(un);
1033 
1034  ctx2.grad_base_value(t2);
1035  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
1036 
1037  ctx2.base_value(t2p);
1038  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
1039 
1040  ctx1.base_value(t1);
1041  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
1042 
1043  ctx1i.base_value(t1i);
1044  vectorize_base_tensor(t1i, tv1i, ndof1i, pf1i->target_dim(), Q);
1045 
1046  ctx3i.base_value(t3i);
1047  vectorize_base_tensor(t3i, tv3i, ndof3i, pf3i->target_dim(), Q);
1048 
1049  ctx3.base_value(t3);
1050  vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
1051 
1052 
1053  for (size_type i = 0; i < ndof1; ++i) // To be optimized
1054  for (size_type j = 0; j < ndof2; ++j)
1055  for (size_type k1 = 0; k1 < Q; ++k1) {
1056  scalar_type b(0), a = coeff *
1057  (tv1p(i, k1) - (i < ndof1i ? tv1i(i, k1) : 0.));
1058  for (size_type k2 = 0; k2 < N; ++k2)
1059  b += a * tv2.as_vector()[j+(k1 + k2*Q)*ndof2] * un[k2];
1060  M1(j, i) += b;
1061  }
1062 
1063  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1064  for (size_type j = 0; j < ndof3; ++j)
1065  for (size_type k = 0; k < Q; ++k)
1066  M7(i, j) += coeff * normun * tv3p(i,k) * tv3p(j, k);
1067 
1068  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1069  for (size_type j = 0; j < ndof2; ++j)
1070  for (size_type k = 0; k < Q; ++k)
1071  M8(i, j) += coeff * normun * tv3p(i,k) * tv2p(j, k);
1072 
1073  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1074  for (size_type j = 0; j < ndof1; ++j)
1075  for (size_type k = 0; k < Q; ++k)
1076  M9(i, j) += coeff * normun * tv3p(i, k) * tv1p(j, k);
1077 
1078  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1079  for (size_type j = 0; j < ndof3i; ++j)
1080  for (size_type k = 0; k < Q; ++k)
1081  M10(i, j) += coeff * normun * tv3p(i,k) * tv3i(j, k);
1082  }
1083  }
1084 
1085  // Add the constraint with penalization
1086  scalar_type coeff_p = pow(area, -1. - 2./scalar_type(N));
1087  gmm::mult(gmm::transposed(M4), M4, aux2);
1088  gmm::add (gmm::scaled(aux2, coeff_p), M2);
1089  gmm::mult(gmm::transposed(M4), M3, aux1);
1090  gmm::add (gmm::scaled(aux1, coeff_p), M1);
1091 
1092  if (pf2->target_dim() == 1 && Q > 1) {
1093  gmm::sub_slice I(0, ndof2/Q, Q);
1094  gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
1095  for (size_type i = 0; i < Q; ++i) {
1096  gmm::sub_slice I2(i, ndof2/Q, Q);
1097  gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
1098  }
1099  } else
1100  { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
1101 
1102  if (pf3->target_dim() == 1 && Q > 1) {
1103  gmm::sub_slice I(0, ndof3/Q, Q);
1104  gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
1105  for (size_type i = 0; i < Q; ++i) {
1106  gmm::sub_slice I2(i, ndof3/Q, Q);
1107  gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
1108  }
1109  } else
1110  { gmm::copy(M7, M7inv); gmm::lu_inverse(M7inv); }
1111 
1112  gmm::mult(M2inv, M1, MD);
1113  gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
1114 
1115  // S = (I - inv(M7)*M10)*inv(M7)*(M9 - M8*MD)
1116  base_matrix MPB(ndof3, ndof3);
1117  gmm::mult(M7inv, M10, MPB);
1118  gmm::copy(gmm::identity_matrix(), M10);
1119  gmm::add(gmm::scaled(MPB, scalar_type(-1)), M10);
1120 
1121  base_matrix MPC(ndof3, ndof1);
1122  gmm::mult(gmm::scaled(M8, scalar_type(-1)), MD, MPC);
1123  gmm::add(M9, MPC);
1124  gmm::mult(M7inv, MPC, M9);
1125  gmm::mult(M10, M9, M);
1126  gmm::clean(M, 1E-13);
1127  }
1128  };
1129 
1130  void add_HHO_stabilization(model &md, std::string name) {
1131  pelementary_transformation
1132  p = std::make_shared<_HHO_stabilization>();
1133  md.add_elementary_transformation(name, p);
1134  }
1135 
1136 #endif
1137 
1138  class _HHO_symmetrized_stabilization
1139  : public virtual_elementary_transformation {
1140 
1141  public:
1142 
1143  virtual void give_transformation(const mesh_fem &mf1, const mesh_fem &mf2,
1144  size_type cv, base_matrix &M) const {
1145  // The reconstructed variable "S" is described on mf2 and computed by
1146  // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)))
1147  // where P_{\dT} et P_T are L2 projections on the boundary and on the
1148  // interior of T on the corresponding discrete spaces.
1149  // Note that P_{\dT}(v_{dT}) = v_{dT} and P_T(v_T) = v_T and D is
1150  // the reconstructed value on P^{k+1} given by the formula:
1151  // \int_T Sym(Grad(D)).Grad(w) = \int_T Sym(Grad(v_T)).Grad(w)
1152  // + \int_{dT}(v_{dT} - v_T).(Sym(Grad(w)).n)
1153  // with the constraints
1154  // \int_T D = \int_T v_T
1155  // \int_T Skew(Grad(D)) = 0.5\int_{dT}(n x v_{dT} - v_{dT} x n)
1156  // where "w" is the test function arbitrary in mf2, "v_T" is the field
1157  // inside the element whose gradient is to be reconstructed,
1158  // "v_{dT}" is the field on the boundary of T and "n" is the outward
1159  // unit normal.
1160  // The implemented formula is
1161  // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)) )
1162  // by the mean of the projection matrix from P^{k+1} to the target space
1163  // and the projection matrix from interior space to the boundary space.
1164  // As it is built, S(v) is zero on interior dofs.
1165 
1166  // Obtaining the fem descriptors
1167  pfem pf1 = mf1.fem_of_element(cv);
1168  short_type degree = pf1->estimated_degree();
1169  bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
1170  pfem pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be changed to an
1171  // interior PK method
1172  pfem pf3 = mf2.fem_of_element(cv);
1173  pfem pf1i = interior_fem_of_hho_method(pf1);
1174  pfem pf3i = interior_fem_of_hho_method(pf3);
1175 
1176  papprox_integration pim
1177  = classical_approx_im(pgt, dim_type(2*degree+2))->approx_method();
1178 
1179  base_matrix G;
1180  bgeot::vectors_to_base_matrix(G, mf1.linked_mesh().points_of_convex(cv));
1181 
1182  bgeot::pgeotrans_precomp pgp
1183  = HHO_pgp_pool(pgt, pim->pintegration_points());
1184  pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
1185  pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
1186  pfem_precomp pfp3 = HHO_pfp_pool(pf3, pim->pintegration_points());
1187  pfem_precomp pfp1i = HHO_pfp_pool(pf1i, pim->pintegration_points());
1188  pfem_precomp pfp3i = HHO_pfp_pool(pf3i, pim->pintegration_points());
1189 
1190  fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
1191  fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
1192  fem_interpolation_context ctx3(pgp, pfp3, 0, G, cv);
1193  fem_interpolation_context ctx1i(pgp, pfp1i, 0, G, cv);
1194  fem_interpolation_context ctx3i(pgp, pfp3i, 0, G, cv);
1195 
1196  size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
1197  GMM_ASSERT1(Q == N, "This transformation works only for vector fields "
1198  "having the same dimension as the domain");
1199  base_vector un(N);
1200  size_type qmult1 = N / pf1->target_dim();
1201  size_type ndof1 = pf1->nb_dof(cv) * qmult1;
1202  size_type qmult2 = N / pf2->target_dim();
1203  size_type ndof2 = pf2->nb_dof(cv) * qmult2;
1204  size_type qmult3 = Q / pf3->target_dim();
1205  size_type ndof3 = pf3->nb_dof(cv) * qmult3;
1206  size_type qmult1i = Q / pf1i->target_dim();
1207  size_type ndof1i = pf1i->nb_dof(cv) * qmult1i;
1208  size_type qmult3i = Q / pf3i->target_dim();
1209  size_type ndof3i = pf3i->nb_dof(cv) * qmult3i;
1210 
1211 
1212  base_tensor t1, t2, t3, t1i, t3i, tv1, tv2, t1p, t2p;
1213  base_matrix tv1p, tv2p, tv3p, tv1i, tv3i;
1214  base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
1215  base_matrix M3(N, ndof1), M4(N, ndof2);
1216  base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
1217  base_matrix M5(N*N, ndof1), M6(N*N, ndof2);
1218  base_matrix M7(ndof3, ndof3), M7inv(ndof3, ndof3), M8(ndof3, ndof2);
1219  base_matrix M9(ndof3, ndof1), M10(ndof3, ndof3), MD(ndof2, ndof1);
1220  scalar_type area(0);
1221 
1222  // Integrals on the element : \int_T Sym(Grad(D)).Grad(w) (M2)
1223  // \int_T Sym(Grad(v_T)).Grad(w) (M1)
1224  // \int_T D (M4) and \int_T v_T (M3)
1225  // \int_T Skew(Grad(D)) (M6)
1226  for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
1227  ctx1.set_ii(ipt); ctx2.set_ii(ipt); ctx3.set_ii(ipt);
1228  scalar_type coeff = pim->coeff(ipt) * ctx1.J();
1229  area += coeff;
1230 
1231  ctx1.grad_base_value(t1);
1232  vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
1233 
1234  ctx1.base_value(t1p);
1235  vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), Q);
1236 
1237  ctx2.grad_base_value(t2);
1238  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
1239 
1240  ctx2.base_value(t2p);
1241  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
1242 
1243  ctx3.base_value(t3);
1244  vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
1245 
1246  for (size_type i = 0; i < ndof2; ++i) // To be optimized
1247  for (size_type j = 0; j < ndof2; ++j)
1248  for (size_type k1 = 0; k1 < N; ++k1)
1249  for (size_type k2 = 0; k2 < N; ++k2)
1250  M2(j, i) += coeff * tv2.as_vector()[i+(k1+k2*N)*ndof2]
1251  * 0.5 * (tv2.as_vector()[j+(k1+k2*N)*ndof2] +
1252  tv2.as_vector()[j+(k2+k1*N)*ndof2]);
1253 
1254  for (size_type i = 0; i < ndof2; ++i)
1255  for (size_type k = 0; k < N; ++k)
1256  M4(k, i) += coeff * tv2p(i, k);
1257 
1258  for (size_type i = 0; i < ndof2; ++i)
1259  for (size_type k1 = 0; k1 < N; ++k1)
1260  for (size_type k2 = 0; k2 < N; ++k2)
1261  M6(k1+k2*N, i) += 0.5*coeff*(tv2.as_vector()[i+(k1+k2*N)*ndof2] -
1262  tv2.as_vector()[i+(k2+k1*N)*ndof2]);
1263 
1264  for (size_type i = 0; i < ndof1; ++i) // To be optimized
1265  for (size_type j = 0; j < ndof2; ++j)
1266  for (size_type k1 = 0; k1 < N; ++k1)
1267  for (size_type k2 = 0; k2 < N; ++k2)
1268  M1(j, i) += coeff * tv1.as_vector()[i+(k1+k2*N)*ndof1]
1269  * 0.5 * (tv2.as_vector()[j+(k1+k2*N)*ndof2] +
1270  tv2.as_vector()[j+(k2+k1*N)*ndof2]);
1271 
1272  for (size_type i = 0; i < ndof1; ++i)
1273  for (size_type k = 0; k < N; ++k)
1274  M3(k, i) += coeff * tv1p(i, k);
1275 
1276  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1277  for (size_type j = 0; j < ndof3; ++j)
1278  for (size_type k = 0; k < N; ++k)
1279  M7(i, j) += coeff * tv3p(i,k) * tv3p(j, k);
1280 
1281  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1282  for (size_type j = 0; j < ndof2; ++j)
1283  for (size_type k = 0; k < N; ++k)
1284  M8(i, j) += coeff * tv3p(i,k) * tv2p(j, k);
1285 
1286  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1287  for (size_type j = 0; j < ndof1; ++j)
1288  for (size_type k = 0; k < Q; ++k)
1289  M9(i, j) += coeff * tv3p(i, k) * tv1p(j, k);
1290  }
1291 
1292  // Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
1293  // \int_{dT} Skew(n x v_{dT} - v_{dT} x n) (M5)
1294  for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
1295  ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctx3.set_face_num(ifc);
1296  ctx1i.set_face_num(ifc); ctx3i.set_face_num(ifc);
1297  size_type first_ind = pim->ind_first_point_on_face(ifc);
1298  for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
1299  ctx1.set_ii(first_ind+ipt); ctx2.set_ii(first_ind+ipt);
1300  ctx3.set_ii(first_ind+ipt); ctx1i.set_ii(first_ind+ipt);
1301  ctx3i.set_ii(first_ind+ipt);
1302  scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
1303  gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
1304  scalar_type normun = gmm::vect_norm2(un);
1305 
1306  ctx2.grad_base_value(t2);
1307  vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
1308 
1309  ctx2.base_value(t2p);
1310  vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
1311 
1312  ctx1.base_value(t1);
1313  vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
1314 
1315  ctx1i.base_value(t1i);
1316  vectorize_base_tensor(t1i, tv1i, ndof1i, pf1i->target_dim(), Q);
1317 
1318  ctx3i.base_value(t3i);
1319  vectorize_base_tensor(t3i, tv3i, ndof3i, pf3i->target_dim(), Q);
1320 
1321  ctx3.base_value(t3);
1322  vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
1323 
1324  for (size_type i = 0; i < ndof1; ++i) // To be optimized
1325  for (size_type j = 0; j < ndof2; ++j)
1326  for (size_type k1 = 0; k1 < N; ++k1) {
1327  scalar_type b(0), a = coeff *
1328  (tv1p(i, k1) - (i < ndof1i ? tv1i(i, k1) : 0.));
1329  for (size_type k2 = 0; k2 < N; ++k2)
1330  b += a * 0.5 * (tv2.as_vector()[j+(k1 + k2*N)*ndof2] +
1331  tv2.as_vector()[j+(k2 + k1*N)*ndof2])* un[k2];
1332  M1(j, i) += b;
1333  }
1334 
1335  for (size_type i = 0; i < ndof1; ++i)
1336  for (size_type k1 = 0; k1 < N; ++k1)
1337  for (size_type k2 = 0; k2 < N; ++k2)
1338  M5(k1+k2*N, i) += 0.5 * coeff * (tv1p(i, k1) * un[k2] -
1339  tv1p(i, k2) * un[k1]);
1340 
1341  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1342  for (size_type j = 0; j < ndof3; ++j)
1343  for (size_type k = 0; k < N; ++k)
1344  M7(i, j) += coeff * normun * tv3p(i,k) * tv3p(j, k);
1345 
1346  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1347  for (size_type j = 0; j < ndof2; ++j)
1348  for (size_type k = 0; k < N; ++k)
1349  M8(i, j) += coeff * normun * tv3p(i,k) * tv2p(j, k);
1350 
1351  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1352  for (size_type j = 0; j < ndof1; ++j)
1353  for (size_type k = 0; k < Q; ++k)
1354  M9(i, j) += coeff * normun * tv3p(i, k) * tv1p(j, k);
1355 
1356  for (size_type i = 0; i < ndof3; ++i) // To be optimized
1357  for (size_type j = 0; j < ndof3i; ++j)
1358  for (size_type k = 0; k < N; ++k)
1359  M10(i, j) += coeff * normun * tv3p(i,k) * tv3i(j, k);
1360  }
1361  }
1362 
1363  // Add the constraint with penalization
1364  scalar_type coeff_p1 = pow(area, -1. - 2./scalar_type(N));
1365  scalar_type coeff_p2 = pow(area, -1. - 1./scalar_type(N));
1366  gmm::mult(gmm::transposed(M4), M4, aux2);
1367  gmm::add (gmm::scaled(aux2, coeff_p1), M2);
1368  gmm::mult(gmm::transposed(M6), M6, aux2);
1369  gmm::add (gmm::scaled(aux2, coeff_p2), M2);
1370  gmm::mult(gmm::transposed(M4), M3, aux1);
1371  gmm::add (gmm::scaled(aux1, coeff_p1), M1);
1372  gmm::mult(gmm::transposed(M6), M5, aux1);
1373  gmm::add (gmm::scaled(aux1, coeff_p2), M1);
1374 
1375  gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv);
1376 
1377  if (pf3->target_dim() == 1 && Q > 1) {
1378  gmm::sub_slice I(0, ndof3/Q, Q);
1379  gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
1380  for (size_type i = 0; i < Q; ++i) {
1381  gmm::sub_slice I2(i, ndof3/Q, Q);
1382  gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
1383  }
1384  } else
1385  { gmm::copy(M7, M7inv); gmm::lu_inverse(M7inv); }
1386 
1387  gmm::mult(M2inv, M1, MD);
1388  gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
1389 
1390  // S = (I - inv(M7)*M10)*inv(M7)*(M9 - M8*MD)
1391  base_matrix MPB(ndof3, ndof3);
1392  gmm::mult(M7inv, M10, MPB);
1393  gmm::copy(gmm::identity_matrix(), M10);
1394  gmm::add(gmm::scaled(MPB, scalar_type(-1)), M10);
1395 
1396  base_matrix MPC(ndof3, ndof1);
1397  gmm::mult(gmm::scaled(M8, scalar_type(-1)), MD, MPC);
1398  gmm::add(M9, MPC);
1399  gmm::mult(M7inv, MPC, M9);
1400  gmm::mult(M10, M9, M);
1401  gmm::clean(M, 1E-13);
1402  }
1403  };
1404 
1405  void add_HHO_symmetrized_stabilization(model &md, std::string name) {
1406  pelementary_transformation
1407  p = std::make_shared<_HHO_symmetrized_stabilization>();
1408  md.add_elementary_transformation(name, p);
1409  }
1410 
1411 
1412 
1413 
1414 
1415 } /* end of namespace getfem. */
1416 
bgeot::geotrans_precomp_pool
The object geotrans_precomp_pool Allow to allocate a certain number of geotrans_precomp and automatic...
Definition: bgeot_geometric_trans.h:382
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:114
getfem::model::add_elementary_transformation
void add_elementary_transformation(const std::string &name, pelementary_transformation ptrans)
Add an elementary transformation to the model to be used with the generic assembly.
Definition: getfem_models.h:1132
getfem_HHO.h
Tools for Hybrid-High-Order methods.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:556
gmm::clean
void clean(L &l, double threshold)
Clean a vector or matrix (replace near-zero entries with zeroes).
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
gmm::lu_inverse
void lu_inverse(const DenseMatrixLU &LU, const Pvector &pvector, const DenseMatrix &AInv_)
Given an LU factored matrix, build the inverse of the matrix.
Definition: gmm_dense_lu.h:211
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
getfem::classical_fem
pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k, bool complete=false)
Give a pointer on the structures describing the classical polynomial fem of degree k on a given conve...
Definition: getfem_fem.cc:4566
getfem::interior_fem_of_hho_method
pfem interior_fem_of_hho_method(pfem hho_method)
Specific function for a HHO method to obtain the method in the interior.
Definition: getfem_fem_composite.cc:865
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:185
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::add_HHO_symmetrized_stabilization
void add_HHO_symmetrized_stabilization(model &md, std::string name)
Add to the model the elementary transformation corresponding to the HHO stabilization operator using ...
Definition: getfem_HHO.cc:1405
getfem::classical_approx_im
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree)
try to find an approximate integration method for the geometric transformation pgt which is able to i...
Definition: getfem_integration.cc:1287
getfem::add_HHO_reconstructed_value
void add_HHO_reconstructed_value(model &md, std::string name)
Add to the model the elementary transformation corresponding to the reconstruction of the variable.
Definition: getfem_HHO.cc:458
getfem::add_HHO_reconstructed_symmetrized_value
void add_HHO_reconstructed_symmetrized_value(model &md, std::string name)
Add to the model the elementary transformation corresponding to the reconstruction of the variable us...
Definition: getfem_HHO.cc:640
getfem::add_HHO_reconstructed_symmetrized_gradient
void add_HHO_reconstructed_symmetrized_gradient(model &md, std::string name)
Add to the model the elementary transformation corresponding to the reconstruction of a symmetrized g...
Definition: getfem_HHO.cc:294
getfem::add_HHO_stabilization
void add_HHO_stabilization(model &md, std::string name)
Add to the model the elementary transformation corresponding to the HHO stabilization operator.
Definition: getfem_HHO.cc:1130
getfem::add_HHO_reconstructed_gradient
void add_HHO_reconstructed_gradient(model &md, std::string name)
Add to the model the elementary transformation corresponding to the reconstruction of a gradient for ...
Definition: getfem_HHO.cc:158