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