GetFEM  5.5
getfem_fem_composite.cc
1 /*===========================================================================
2 
3  Copyright (C) 2002-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/bgeot_poly_composite.h"
23 #include "getfem/getfem_integration.h"
24 #include "getfem/getfem_mesh_fem.h"
25 #include "getfem/dal_naming_system.h"
26 
27 namespace getfem {
28 
29  typedef const fem<bgeot::polynomial_composite> * ppolycompfem;
30 
31  struct polynomial_composite_fem : public fem<bgeot::polynomial_composite> {
32  mesh m;
33  mutable bgeot::mesh_precomposite mp;
34  mesh_fem mf;
35  mutable bgeot::pgeotrans_precomp pgp;
36  mutable bgeot::pgeometric_trans pgt_stored;
37  bgeot::pstored_point_tab pspt;
38 
39 
40  virtual void mat_trans(base_matrix &M, const base_matrix &G,
41  bgeot::pgeometric_trans pgt) const;
42 
43 
44  polynomial_composite_fem(const mesh &m_, const mesh_fem &mf_)
45  : m(m_), mp(m), mf(m), pgp(0), pgt_stored(0) {
46  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv)
47  mf.set_finite_element(cv, mf_.fem_of_element(cv));
48  mf.nb_dof();
49  pspt = store_point_tab(m.points());
50  // verification for the non-equivalent fems
51  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
52  pfem pf1 = mf.fem_of_element(cv);
53  if (!(pf1->is_equivalent())) {
54  dal::bit_vector unshareable;
55  for (const auto &i : mf.ind_basic_dof_of_element(cv))
56  unshareable.add(i);
57 
58  for (dal::bv_visitor cv2(mf.convex_index()); !cv2.finished(); ++cv2) {
59  if (cv2 != cv)
60  for (const auto &i : mf.ind_basic_dof_of_element(cv2))
61  GMM_ASSERT1(!(unshareable.is_in(i)), "Non equivalent elements "
62  "are allowed only if they do not share their dofs");
63  }
64  }
65  }
66  }
67  };
68 
69  void polynomial_composite_fem::mat_trans(base_matrix &M, const base_matrix &G,
70  bgeot::pgeometric_trans pgt) const {
71  dim_type N = dim_type(G.nrows());
72  gmm::copy(gmm::identity_matrix(), M);
73  base_matrix G1, M1;
74 
75  if (pgt != pgt_stored) {
76  pgt_stored = pgt;
77  pgp = bgeot::geotrans_precomp(pgt, pspt, 0);
78  }
79 
80  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
81  pfem pf1 = mf.fem_of_element(cv);
82  if (!(pf1->is_equivalent())) {
83  size_type npt=m.nb_points_of_convex(cv);
84  size_type ndof=mf.nb_basic_dof_of_element(cv);
85  GMM_ASSERT1(ndof == pf1->nb_base(0) && ndof == pf1->nb_dof(0),
86  "Sorry, limited implementation for the moment");
87  // Compute the local G
88  G1.resize(npt, N);
89  M1.resize(ndof, ndof);
90  for (size_type i = 0; i < npt; ++i)
91  gmm::copy(pgp->transform(m.ind_points_of_convex(cv)[i], G),
92  gmm::mat_col(G1, i));
93  // Call for the local M
94  pf1->mat_trans(M1, G1, m.trans_of_convex(cv));
95  gmm::sub_index I(mf.ind_basic_dof_of_element(cv));
96  // I is in fact an interval ...
97  gmm::copy(M1, gmm::sub_matrix(M, I, I));
98  }
99  }
100  }
101 
102  static pfem composite_fe_method(const getfem::mesh &m,
103  const mesh_fem &mf, bgeot::pconvex_ref cr,
104  bool ff = false) {
105 
106  GMM_ASSERT1(!mf.is_reduced(),
107  "Sorry, does not work for reduced mesh_fems");
108  auto p = std::make_shared<polynomial_composite_fem>(m, mf);
109 
110  p->mref_convex() = cr;
111  p->dim() = cr->structure()->dim();
112  p->is_polynomialcomp() = p->is_equivalent() = p->is_standard() = true;
113  p->is_polynomial() = false;
114  p->is_lagrange() = true;
115  p->estimated_degree() = 0;
116  p->init_cvs_node();
117 
118  std::vector<bgeot::polynomial_composite> base(mf.nb_basic_dof());
119  std::fill(base.begin(), base.end(),
120  bgeot::polynomial_composite(p->mp, true, ff));
121  std::vector<pdof_description> dofd(mf.nb_basic_dof());
122 
123  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
124  pfem pf1 = mf.fem_of_element(cv);
125  if (!pf1->is_lagrange()) p->is_lagrange() = false;
126  if (!(pf1->is_equivalent())) p->is_equivalent() = false;
127 
128  GMM_ASSERT1(pf1->is_polynomial(), "Only for polynomial fems.");
129  ppolyfem pf = ppolyfem(pf1.get());
130  p->estimated_degree() = std::max(p->estimated_degree(),
131  pf->estimated_degree());
132  for (size_type k = 0; k < pf->nb_dof(cv); ++k) {
133  size_type igl = mf.ind_basic_dof_of_element(cv)[k];
134  base_poly fu = pf->base()[k];
135  base[igl].set_poly_of_subelt(cv, fu);
136  dofd[igl] = pf->dof_types()[k];
137  }
138  }
139  p->base().resize(mf.nb_basic_dof());
140  for (size_type k = 0; k < mf.nb_basic_dof(); ++k) {
141  p->add_node(dofd[k], mf.point_of_basic_dof(k));
142  p->base()[k] = base[k];
143  }
144  return p;
145  }
146 
147  typedef dal::naming_system<virtual_fem>::param_list fem_param_list;
148 
149  pfem structured_composite_fem_method(fem_param_list &params,
150  std::vector<dal::pstatic_stored_object> &dependencies) {
151  GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
152  << params.size() << " should be 2.");
153  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 0,
154  "Bad type of parameters");
155  pfem pf = params[0].method();
156  int k = int(::floor(params[1].num() + 0.01));
157  GMM_ASSERT1(((pf->is_polynomial()) || !(pf->is_equivalent())) && k > 0
158  && k <= 150 && double(k) == params[1].num(), "Bad parameters");
159  bgeot::pbasic_mesh pm;
160  bgeot::pmesh_precomposite pmp;
161 
162  structured_mesh_for_convex(pf->ref_convex(0), short_type(k), pm, pmp);
163 
164  mesh m(*pm);
165  mesh_fem mf(m);
166  mf.set_finite_element(pm->convex_index(), pf);
167  pfem p = composite_fe_method(m, mf, pf->ref_convex(0));
168  dependencies.push_back(p->ref_convex(0));
169  dependencies.push_back(p->node_tab(0));
170  return p;
171  }
172 
173  pfem PK_composite_hierarch_fem(fem_param_list &params,
174  std::vector<dal::pstatic_stored_object> &) {
175  GMM_ASSERT1(params.size() == 3, "Bad number of parameters : "
176  << params.size() << " should be 3.");
177  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
178  params[2].type() == 0, "Bad type of parameters");
179  int n = int(::floor(params[0].num() + 0.01));
180  int k = int(::floor(params[1].num() + 0.01));
181  int s = int(::floor(params[2].num() + 0.01)), t;
182  GMM_ASSERT1(n > 0 && n < 100 && k > 0 && k <= 150 && s > 0 && s <= 150 &&
183  (!(s & 1) || (s == 1)) && double(s) == params[2].num() &&
184  double(n) == params[0].num() && double(k) == params[1].num(),
185  "Bad parameters");
186  std::stringstream name;
187  if (s == 1)
188  name << "FEM_STRUCTURED_COMPOSITE(FEM_PK(" << n << "," << k << "),1)";
189  else {
190  for (t = 2; t <= s; ++t) if ((s % t) == 0) break;
191  name << "FEM_GEN_HIERARCHICAL(FEM_PK_HIERARCHICAL_COMPOSITE(" << n
192  << "," << k << "," << s/t << "), FEM_STRUCTURED_COMPOSITE(FEM_PK("
193  << n << "," << k << ")," << s << "))";
194  }
195  return fem_descriptor(name.str());
196  }
197 
198  pfem PK_composite_full_hierarch_fem(fem_param_list &params,
199  std::vector<dal::pstatic_stored_object> &) {
200  GMM_ASSERT1(params.size() == 3, "Bad number of parameters : "
201  << params.size() << " should be 3.");
202  GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
203  params[2].type() == 0, "Bad type of parameters");
204  int n = int(::floor(params[0].num() + 0.01));
205  int k = int(::floor(params[1].num() + 0.01));
206  int s = int(::floor(params[2].num() + 0.01)), t;
207  GMM_ASSERT1(n > 0 && n < 100 && k > 0 && k <= 150 && s > 0 && s <= 150 &&
208  (!(s & 1) || (s == 1)) && double(s) == params[2].num() &&
209  double(n) == params[0].num() && double(k) == params[1].num(),
210  "Bad parameters");
211  std::stringstream name;
212  if (s == 1)
213  name << "FEM_STRUCTURED_COMPOSITE(FEM_PK_HIERARCHICAL(" << n << ","
214  << k << "),1)";
215  else {
216  for (t = 2; t <= s; ++t) if ((s % t) == 0) break;
217  name << "FEM_GEN_HIERARCHICAL(FEM_PK_FULL_HIERARCHICAL_COMPOSITE(" << n
218  << "," << k << "," << s/t
219  << "), FEM_STRUCTURED_COMPOSITE(FEM_PK_HIERARCHICAL("
220  << n << "," << k << ")," << s << "))";
221  }
222  return fem_descriptor(name.str());
223  }
224 
225 
226  /* ******************************************************************** */
227  /* P1 with piecewise linear bubble function on a triangle. */
228  /* ******************************************************************** */
229 
230  struct P1bubbletriangle__ : public fem<bgeot::polynomial_composite> {
231  mesh m;
232  bgeot::mesh_precomposite mp;
233  P1bubbletriangle__(void);
234  };
235 
236  P1bubbletriangle__::P1bubbletriangle__(void) {
237 
238  m.clear();
239  size_type i0 = m.add_point(base_node(1.0/3.0, 1.0/3.0));
240  size_type i1 = m.add_point(base_node(0.0, 0.0));
241  size_type i2 = m.add_point(base_node(1.0, 0.0));
242  size_type i3 = m.add_point(base_node(0.0, 1.0));
243  m.add_triangle(i0, i2, i3);
244  m.add_triangle(i0, i3, i1);
245  m.add_triangle(i0, i1, i2);
246  mp.initialise(m);
247 
248  std::stringstream s("1-x-y;1-x-y;1-x-y;x;x;x;y;y;y;3-3*x-3*y;3*x;3*y;");
249 
250  bgeot::pconvex_ref cr = bgeot::simplex_of_reference(2);
251  mref_convex() = cr;
252  dim() = cr->structure()->dim();
253  is_polynomialcomp() = true;
254  is_equivalent() = true;
255  is_polynomial() = false;
256  is_lagrange() = false;
257  is_standard() = true;
258  estimated_degree() = 3;
259  init_cvs_node();
260 
261  base()=std::vector<bgeot::polynomial_composite>
262  (4, bgeot::polynomial_composite(mp, false));
263  for (size_type k = 0; k < 4; ++k)
264  for (size_type ic = 0; ic < 3; ++ic)
265  base()[k].set_poly_of_subelt(ic, bgeot::read_base_poly(2, s));
266 
267  for (size_type i = 0; i < 3; ++i) {
268  base_node pt(0.0, 0.0);
269  if (i) pt[i-1] = 1.0;
270  add_node(lagrange_dof(2), pt);
271  }
272 
273  add_node(bubble1_dof(2), base_node(1.0/3.0, 1.0/3.0));
274  }
275 
276 
277  pfem P1bubbletriangle_fem
278  (fem_param_list &params,
279  std::vector<dal::pstatic_stored_object> &dependencies) {
280  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
281  << params.size() << " should be 0.");
282  pfem p = std::make_shared<P1bubbletriangle__>();
283  dependencies.push_back(p->ref_convex(0));
284  dependencies.push_back(p->node_tab(0));
285  return p;
286  }
287 
288  /* ******************************************************************** */
289  /* Hsieh-Clough-Tocher C^1 element (composite P3) */
290  /* ******************************************************************** */
291 
292  struct HCT_triangle__ : public fem<bgeot::polynomial_composite> {
293  virtual void mat_trans(base_matrix &M, const base_matrix &G,
294  bgeot::pgeometric_trans pgt) const;
295  mesh m;
296  mutable bgeot::base_small_vector true_normals[3];
297  mutable bgeot::mesh_precomposite mp;
298  mutable bgeot::pgeotrans_precomp pgp;
299  mutable pfem_precomp pfp;
300  mutable bgeot::pgeometric_trans pgt_stored;
301  mutable base_matrix K;
302 
303  HCT_triangle__(void);
304  };
305 
306  void HCT_triangle__::mat_trans(base_matrix &M, const base_matrix &G,
307  bgeot::pgeometric_trans pgt) const {
308 
309  dim_type N = dim_type(G.nrows());
310 
311  GMM_ASSERT1(N == 2, "Sorry, this version of HCT "
312  "element works only on dimension two.");
313  if (pgt != pgt_stored) {
314  pgt_stored = pgt;
315  pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
316  pfp = fem_precomp(std::make_shared<HCT_triangle__>(), node_tab(0), 0);
317  }
318  gmm::copy(gmm::identity_matrix(), M);
319 
320  gmm::mult(G, pgp->grad(0), K);
321  for (size_type i = 0; i < 3; ++i) {
322  if (i && !(pgt->is_linear())) gmm::mult(G, pgp->grad(3*i), K);
323  gmm::copy(K, gmm::sub_matrix(M, gmm::sub_interval(1+3*i, 2)));
324  }
325 
326  // take the normal derivatives into account
327  static base_matrix W(3, 12);
328  base_small_vector norient(M_PI, M_PI * M_PI);
329  if (pgt->is_linear()) gmm::lu_inverse(K);
330  for (unsigned i = 9; i < 12; ++i) {
331  if (!(pgt->is_linear()))
332  { gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
333  bgeot::base_small_vector n(2), v(2);
334  gmm::mult(gmm::transposed(K), cvr->normals()[i-9], n);
335  n /= gmm::vect_norm2(n);
336 
337  scalar_type ps = gmm::vect_sp(n, norient);
338  if (ps < 0) n *= scalar_type(-1);
339  true_normals[i-9] = n;
340 
341  if (gmm::abs(ps) < 1E-8)
342  GMM_WARNING2("HCT_triangle : "
343  "The normal orientation may be not correct");
344  gmm::mult(K, n, v);
345  const bgeot::base_tensor &t = pfp->grad(i);
346  // cout << "t = " << t << endl;
347  for (unsigned j = 0; j < 12; ++j)
348  W(i-9, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1];
349  }
350 
351  static base_matrix A(3, 3);
352  static bgeot::base_vector w(3), coeff(3);
353  static gmm::sub_interval SUBI(9, 3), SUBJ(0, 3);
354  gmm::copy(gmm::sub_matrix(W, SUBJ, SUBI), A);
355  gmm::lu_inverse(A);
356  gmm::copy(gmm::transposed(A), gmm::sub_matrix(M, SUBI));
357 
358  for (unsigned j = 0; j < 9; ++j) {
359  gmm::mult(W, gmm::mat_row(M, j), w);
360  gmm::mult(A, gmm::scaled(w, -1.0), coeff);
361  gmm::copy(coeff, gmm::sub_vector(gmm::mat_row(M, j), SUBI));
362  }
363  }
364 
365  HCT_triangle__::HCT_triangle__(void) : pgt_stored(0), K(2, 2) {
366 
367  m.clear();
368  size_type i0 = m.add_point(base_node(1.0/3.0, 1.0/3.0));
369  size_type i1 = m.add_point(base_node(0.0, 0.0));
370  size_type i2 = m.add_point(base_node(1.0, 0.0));
371  size_type i3 = m.add_point(base_node(0.0, 1.0));
372  m.add_triangle(i0, i2, i3);
373  m.add_triangle(i0, i3, i1);
374  m.add_triangle(i0, i1, i2);
375  mp.initialise(m);
376 
377 
378  std::stringstream s
379  ("-1 + 9*x + 9*y - 15*x^2 - 30*x*y - 15*y^2 + 7*x^3 + 21*x^2*y + 21*x*y^2 + 7*y^3;"
380  "1 - 3*x^2 - 3*y^2 + 3*x^3 - 3*x^2*y + 2*y^3;"
381  "1 - 3*x^2 - 3*y^2 + 2*x^3 - 3*x*y^2 + 3*y^3;"
382  "-1/6 + 5/2*x - 9/2*x^2 - 4*x*y + 1/2*y^2 + 13/6*x^3 + 4*x^2*y + 3/2*x*y^2 - 1/3*y^3;"
383  "x - 1/2*x^2 - 3*x*y - 7/6*x^3 + 2*x^2*y + 2*x*y^2;"
384  "x - 2*x^2 - 3/2*y^2 + x^3 - 1/2*x*y^2 + 7/3*y^3;"
385  "-1/6 + 5/2*y + 1/2*x^2 - 4*x*y - 9/2*y^2 - 1/3*x^3 + 3/2*x^2*y + 4*x*y^2 + 13/6*y^3;"
386  "y - 3/2*x^2 - 2*y^2 + 7/3*x^3 - 1/2*x^2*y + y^3;"
387  "y - 3*x*y - 1/2*y^2 + 2*x^2*y + 2*x*y^2 - 7/6*y^3;"
388  "1 - 9/2*x - 9/2*y + 9*x^2 + 15*x*y + 6*y^2 - 9/2*x^3 - 21/2*x^2*y - 21/2*x*y^2 - 5/2*y^3;"
389  "3*x^2 - 5/2*x^3 + 3/2*x^2*y;"
390  "3*x^2 - 2*x^3 + 3/2*x*y^2 - 1/2*y^3;"
391  "-1/6 + 3/4*x + 3/4*y - 2*x^2 - 5/2*x*y - y^2 + 17/12*x^3 + 7/4*x^2*y + 7/4*x*y^2 + 5/12*y^3;"
392  "-x^2 + 13/12*x^3 - 1/4*x^2*y;"
393  "-x^2 + x^3 - 1/4*x*y^2 + 1/12*y^3;"
394  "2/3 - 13/4*x - 11/4*y + 9/2*x^2 + 19/2*x*y + 7/2*y^2 - 23/12*x^3 - 23/4*x^2*y - 25/4*x*y^2 - 17/12*y^3;"
395  "-1/2*x^2 + 5/12*x^3 + 9/4*x^2*y;"
396  "-x*y + 1/2*y^2 + 2*x^2*y + 7/4*x*y^2 - 13/12*y^3;"
397  "1 - 9/2*x - 9/2*y + 6*x^2 + 15*x*y + 9*y^2 - 5/2*x^3 - 21/2*x^2*y - 21/2*x*y^2 - 9/2*y^3;"
398  "3*y^2 - 1/2*x^3 + 3/2*x^2*y - 2*y^3;"
399  "3*y^2 + 3/2*x*y^2 - 5/2*y^3;"
400  "2/3 - 11/4*x - 13/4*y + 7/2*x^2 + 19/2*x*y + 9/2*y^2 - 17/12*x^3 - 25/4*x^2*y - 23/4*x*y^2 - 23/12*y^3;"
401  "1/2*x^2 - x*y - 13/12*x^3 + 7/4*x^2*y + 2*x*y^2;"
402  "-1/2*y^2 + 9/4*x*y^2 + 5/12*y^3;"
403  "-1/6 + 3/4*x + 3/4*y - x^2 - 5/2*x*y - 2*y^2 + 5/12*x^3 + 7/4*x^2*y + 7/4*x*y^2 + 17/12*y^3;"
404  "-y^2 + 1/12*x^3 - 1/4*x^2*y + y^3;"
405  "-y^2 - 1/4*x*y^2 + 13/12*y^3;"
406  "-sqrt(2)*2/3 + sqrt(2)*3*x + sqrt(2)*3*y - sqrt(2)*4*x^2 - sqrt(2)*10*x*y - sqrt(2)*4*y^2 + sqrt(2)*5/3*x^3 + sqrt(2)*7*x^2*y + sqrt(2)*7*x*y^2 + sqrt(2)*5/3*y^3;"
407  "sqrt(2)*1/3*x^3 - sqrt(2)*x^2*y;"
408  "-sqrt(2)*x*y^2 + sqrt(2)*1/3*y^3;"
409  "2/3 - 2*x - 4*y + 2*x^2 + 8*x*y + 6*y^2 - 2/3*x^3 - 4*x^2*y - 6*x*y^2 - 8/3*y^3;"
410  "2*x^2 - 4*x*y - 10/3*x^3 + 4*x^2*y + 4*x*y^2;"
411  "-2*y^2 + 2*x*y^2 + 8/3*y^3;"
412  "2/3 - 4*x - 2*y + 6*x^2 + 8*x*y + 2*y^2 - 8/3*x^3 - 6*x^2*y - 4*x*y^2 - 2/3*y^3;"
413  "-2*x^2 + 8/3*x^3 + 2*x^2*y;"
414  "-4*x*y + 2*y^2 + 4*x^2*y + 4*x*y^2 - 10/3*y^3;");
415 
416  bgeot::pconvex_ref cr = bgeot::simplex_of_reference(2);
417  mref_convex() = cr;
418  dim() = cr->structure()->dim();
419  is_polynomialcomp() = true;
420  is_equivalent() = false;
421  is_polynomial() = false;
422  is_lagrange() = false;
423  is_standard() = false;
424  estimated_degree() = 5;
425  init_cvs_node();
426 
427  base()=std::vector<bgeot::polynomial_composite>
428  (12, bgeot::polynomial_composite(mp, false));
429  for (size_type k = 0; k < 12; ++k)
430  for (size_type ic = 0; ic < 3; ++ic)
431  base()[k].set_poly_of_subelt(ic, bgeot::read_base_poly(2, s));
432 
433  for (size_type i = 0; i < 3; ++i) {
434  base_node pt(0.0, 0.0);
435  if (i) pt[i-1] = 1.0;
436  add_node(lagrange_dof(2), pt);
437  add_node(derivative_dof(2, 0), pt);
438  add_node(derivative_dof(2, 1), pt);
439  }
440 
441  add_node(normal_derivative_dof(2), base_node(0.5, 0.5));
442  add_node(normal_derivative_dof(2), base_node(0.0, 0.5));
443  add_node(normal_derivative_dof(2), base_node(0.5, 0.0));
444  }
445 
446  pfem HCT_triangle_fem
447  (fem_param_list &params,
448  std::vector<dal::pstatic_stored_object> &dependencies) {
449  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
450  << params.size() << " should be 0.");
451  pfem p = std::make_shared<HCT_triangle__>();
452  dependencies.push_back(p->ref_convex(0));
453  dependencies.push_back(p->node_tab(0));
454  return p;
455  }
456 
457 
458  /* ******************************************************************** */
459  /* Reduced Hsieh-Clough-Tocher C^1 element (composite P3) */
460  /* ******************************************************************** */
461 
462  struct reduced_HCT_triangle__ : public fem<bgeot::polynomial_composite> {
463  const HCT_triangle__ *HCT;
464  virtual void mat_trans(base_matrix &M, const base_matrix &G,
465  bgeot::pgeometric_trans pgt) const;
466  virtual size_type nb_base(size_type) const { return 12; }
467  mutable base_matrix P, Mhct;
468  reduced_HCT_triangle__(void);
469  };
470 
471  void reduced_HCT_triangle__::mat_trans(base_matrix &M, const base_matrix &G,
472  bgeot::pgeometric_trans pgt) const {
473  HCT->mat_trans(Mhct, G, pgt);
474 
475  P(10, 1)=HCT->true_normals[1][0]*0.5; P(11, 1)=HCT->true_normals[2][0]*0.5;
476  P(10, 2)=HCT->true_normals[1][1]*0.5; P(11, 2)=HCT->true_normals[2][1]*0.5;
477 
478  P( 9, 4)=HCT->true_normals[0][0]*0.5; P(11, 4)=HCT->true_normals[2][0]*0.5;
479  P( 9, 5)=HCT->true_normals[0][1]*0.5; P(11, 5)=HCT->true_normals[2][1]*0.5;
480 
481  P( 9, 7)=HCT->true_normals[0][0]*0.5; P(10, 7)=HCT->true_normals[1][0]*0.5;
482  P( 9, 8)=HCT->true_normals[0][1]*0.5; P(10, 8)=HCT->true_normals[1][1]*0.5;
483 
484  gmm::mult(gmm::transposed(P), Mhct, M);
485  }
486 
487  reduced_HCT_triangle__::reduced_HCT_triangle__(void)
488  : P(12, 9), Mhct(12, 12) {
489  HCT = dynamic_cast<const HCT_triangle__ *>
490  (&(*fem_descriptor("FEM_HCT_TRIANGLE")));
491 
492  bgeot::pconvex_ref cr = bgeot::simplex_of_reference(2);
493  mref_convex() = cr;
494  dim() = cr->structure()->dim();
495  is_polynomialcomp() = true;
496  is_equivalent() = false;
497  is_polynomial() = false;
498  is_lagrange() = false;
499  is_standard() = false;
500  estimated_degree() = 5;
501  base() = HCT->base();
502 
503  gmm::copy(gmm::identity_matrix(), P);
504  init_cvs_node();
505 
506  for (size_type i = 0; i < 3; ++i) {
507  base_node pt(0.0, 0.0);
508  if (i) pt[i-1] = 1.0;
509  add_node(lagrange_dof(2), pt);
510  add_node(derivative_dof(2, 0), pt);
511  add_node(derivative_dof(2, 1), pt);
512  }
513  }
514 
515 
516  pfem reduced_HCT_triangle_fem
517  (fem_param_list &params,
518  std::vector<dal::pstatic_stored_object> &dependencies) {
519  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
520  << params.size() << " should be 0.");
521  pfem p = std::make_shared<reduced_HCT_triangle__>();
522  dependencies.push_back(p->ref_convex(0));
523  dependencies.push_back(p->node_tab(0));
524  return p;
525  }
526 
527 
528  /* ******************************************************************** */
529  /* C1 composite element on quadrilateral (piecewise P3, FVS element). */
530  /* ******************************************************************** */
531 
532  struct quadc1p3__ : public fem<bgeot::polynomial_composite> {
533  virtual void mat_trans(base_matrix &M, const base_matrix &G,
534  bgeot::pgeometric_trans pgt) const;
535  mesh m;
536  mutable bgeot::mesh_precomposite mp;
537  mutable bgeot::pgeotrans_precomp pgp;
538  mutable pfem_precomp pfp;
539  mutable bgeot::pgeometric_trans pgt_stored;
540  mutable base_matrix K;
541  mutable bgeot::base_small_vector true_normals[4];
542  quadc1p3__(void);
543  };
544 
545  void quadc1p3__::mat_trans(base_matrix &M, const base_matrix &G,
546  bgeot::pgeometric_trans pgt) const {
547 
548  dim_type N = dim_type(G.nrows());
549 
550  GMM_ASSERT1(N == 2, "Sorry, this version of reduced HCT "
551  "element works only on dimension two.");
552  if (pgt != pgt_stored) {
553  pgt_stored = pgt;
554  pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
555  pfp = fem_precomp(std::make_shared<quadc1p3__>(), node_tab(0), 0);
556  }
557  gmm::copy(gmm::identity_matrix(), M);
558 
559  gmm::mult(G, pgp->grad(0), K);
560  for (size_type i = 0; i < 4; ++i) {
561  if (i && !(pgt->is_linear())) gmm::mult(G, pgp->grad(3*i), K);
562  gmm::copy(K, gmm::sub_matrix(M, gmm::sub_interval(1+3*i, 2)));
563  }
564 
565  // take the normal derivatives into account
566  static base_matrix W(4, 16);
567  base_small_vector norient(M_PI, M_PI * M_PI);
568  if (pgt->is_linear()) gmm::lu_inverse(K);
569  for (unsigned i = 12; i < 16; ++i) {
570  if (!(pgt->is_linear()))
571  { gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
572  bgeot::base_small_vector n(2), v(2);
573  gmm::mult(gmm::transposed(K), cvr->normals()[i-12], n);
574  n /= gmm::vect_norm2(n);
575 
576  scalar_type ps = gmm::vect_sp(n, norient);
577  if (ps < 0) n *= scalar_type(-1);
578  true_normals[i-12] = n;
579  if (gmm::abs(ps) < 1E-8)
580  GMM_WARNING2("FVS_quadrilateral : "
581  "The normal orientation may be not correct");
582  gmm::mult(K, n, v);
583  const bgeot::base_tensor &t = pfp->grad(i);
584  for (unsigned j = 0; j < 16; ++j)
585  W(i-12, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1];
586  }
587 
588  static base_matrix A(4, 4);
589  static bgeot::base_vector w(4), coeff(4);
590  static gmm::sub_interval SUBI(12, 4), SUBJ(0, 4);
591  gmm::copy(gmm::sub_matrix(W, SUBJ, SUBI), A);
592  gmm::lu_inverse(A);
593  gmm::copy(gmm::transposed(A), gmm::sub_matrix(M, SUBI));
594 
595  for (unsigned j = 0; j < 12; ++j) {
596  gmm::mult(W, gmm::mat_row(M, j), w);
597  gmm::mult(A, gmm::scaled(w, -1.0), coeff);
598  gmm::copy(coeff, gmm::sub_vector(gmm::mat_row(M, j), SUBI));
599  }
600  }
601 
602  quadc1p3__::quadc1p3__(void) : pgt_stored(0), K(2, 2) {
603 
604  m.clear();
605  size_type i0 = m.add_point(base_node(0.0, 0.0));
606  size_type i1 = m.add_point(base_node(1.0, 0.0));
607  size_type i2 = m.add_point(base_node(0.0, 1.0));
608  size_type i3 = m.add_point(base_node(1.0, 1.0));
609  size_type i4 = m.add_point(base_node(0.5, 0.5));
610  m.add_triangle(i1, i3, i4);
611  m.add_triangle(i2, i0, i4);
612  m.add_triangle(i3, i2, i4);
613  m.add_triangle(i0, i1, i4);
614  mp.initialise(m);
615 
616  std::stringstream s
617  ("2 - 3*x - 3*y + 6*x*y + x^3 - 3*x^2*y;"
618  "1 - 3*x^2 - 3*y^2 + x^3 + 3*x^2*y + 2*y^3;"
619  "2 - 3*x - 3*y + 6*x*y - 3*x*y^2 + y^3;"
620  "1 - 3*x^2 - 3*y^2 + 2*x^3 + 3*x*y^2 + y^3;"
621  "1/6 + 1/2*x - y - 3/2*x^2 + 2*x*y + 5/6*x^3 - x^2*y;"
622  "x - 1/2*x^2 - 3*x*y + 1/6*x^3 + x^2*y + 2*x*y^2;"
623  "1/6 + 1/2*x - y - x*y + 3/2*y^2 + 1/2*x*y^2 - 2/3*y^3;"
624  "x - 2*x^2 - 3/2*y^2 + x^3 + 3/2*x*y^2 + 2/3*y^3;"
625  "1/6 - x + 1/2*y + 3/2*x^2 - x*y - 2/3*x^3 + 1/2*x^2*y;"
626  "y - 3/2*x^2 - 2*y^2 + 2/3*x^3 + 3/2*x^2*y + y^3;"
627  "1/6 - x + 1/2*y + 2*x*y - 3/2*y^2 - x*y^2 + 5/6*y^3;"
628  "y - 3*x*y - 1/2*y^2 + 2*x^2*y + x*y^2 + 1/6*y^3;"
629  "-1 + 3*x + 3*y - 6*x*y - 3*y^2 - x^3 + 3*x^2*y + 2*y^3;"
630  "3*x^2 - x^3 - 3*x^2*y;"
631  "-1 + 3*x + 3*y - 6*x*y - 3*y^2 + 3*x*y^2 + y^3;"
632  "3*x^2 - 2*x^3 - 3*x*y^2 + y^3;"
633  "-2/3 + 1/2*x + 2*y - x*y - 2*y^2 + 1/6*x^3 - x^2*y + 2*x*y^2;"
634  "-x^2 + 5/6*x^3 + x^2*y;"
635  "-2/3 + 1/2*x + 2*y - x*y - 2*y^2 + 1/2*x*y^2 + 2/3*y^3;"
636  "-x^2 + x^3 + 3/2*x*y^2 - 2/3*y^3;"
637  "-5/6 + x + 5/2*y + 1/2*x^2 - 3*x*y - 2*y^2 - 2/3*x^3 + 3/2*x^2*y + y^3;"
638  "-1/2*x^2 + 2/3*x^3 + 1/2*x^2*y;"
639  "-5/6 + x + 5/2*y - 2*x*y - 5/2*y^2 + x*y^2 + 5/6*y^3;"
640  "-x*y + 1/2*y^2 + 2*x^2*y - x*y^2 + 1/6*y^3;"
641  "-1 + 3*x + 3*y - 3*x^2 - 6*x*y + x^3 + 3*x^2*y;"
642  "3*y^2 + x^3 - 3*x^2*y - 2*y^3;"
643  "-1 + 3*x + 3*y - 3*x^2 - 6*x*y + 2*x^3 + 3*x*y^2 - y^3;"
644  "3*y^2 - 3*x*y^2 - y^3;"
645  "-5/6 + 5/2*x + y - 5/2*x^2 - 2*x*y + 5/6*x^3 + x^2*y;"
646  "1/2*x^2 - x*y + 1/6*x^3 - x^2*y + 2*x*y^2;"
647  "-5/6 + 5/2*x + y - 2*x^2 - 3*x*y + 1/2*y^2 + x^3 + 3/2*x*y^2 - 2/3*y^3;"
648  "-1/2*y^2 + 1/2*x*y^2 + 2/3*y^3;"
649  "-2/3 + 2*x + 1/2*y - 2*x^2 - x*y + 2/3*x^3 + 1/2*x^2*y;"
650  "-y^2 - 2/3*x^3 + 3/2*x^2*y + y^3;"
651  "-2/3 + 2*x + 1/2*y - 2*x^2 - x*y + 2*x^2*y - x*y^2 + 1/6*y^3;"
652  "-y^2 + x*y^2 + 5/6*y^3;"
653  "1 - 3*x - 3*y + 3*x^2 + 6*x*y + 3*y^2 - x^3 - 3*x^2*y - 2*y^3;"
654  "-x^3 + 3*x^2*y;"
655  "1 - 3*x - 3*y + 3*x^2 + 6*x*y + 3*y^2 - 2*x^3 - 3*x*y^2 - y^3;"
656  "3*x*y^2 - y^3;"
657  "-2/3 + 3/2*x + 2*y - x^2 - 3*x*y - 2*y^2 + 1/6*x^3 + x^2*y + 2*x*y^2;"
658  "5/6*x^3 - x^2*y;"
659  "-2/3 + 3/2*x + 2*y - x^2 - 3*x*y - 2*y^2 + x^3 + 3/2*x*y^2 + 2/3*y^3;"
660  "1/2*x*y^2 - 2/3*y^3;"
661  "-2/3 + 2*x + 3/2*y - 2*x^2 - 3*x*y - y^2 + 2/3*x^3 + 3/2*x^2*y + y^3;"
662  "-2/3*x^3 + 1/2*x^2*y;"
663  "-2/3 + 2*x + 3/2*y - 2*x^2 - 3*x*y - y^2 + 2*x^2*y + x*y^2 + 1/6*y^3;"
664  "-x*y^2 + 5/6*y^3;"
665  "4/3 - 2*x - 4*y + 4*x*y + 4*y^2 + 2/3*x^3 - 4*x*y^2;"
666  "-2/3*x^3;"
667  "4/3 - 2*x - 4*y + 4*x*y + 4*y^2 - 2*x*y^2 - 4/3*y^3;"
668  "-2*x*y^2 + 4/3*y^3;"
669  "-2/3 + 2*x - 2*x^2 + 2/3*x^3;"
670  "2*x^2 - 4*x*y - 2/3*x^3 + 4*x*y^2;"
671  "-2/3 + 2*x - 4*x*y + 2*y^2 + 2*x*y^2 - 4/3*y^3;"
672  "-2*y^2 + 2*x*y^2 + 4/3*y^3;"
673  "4/3 - 4*x - 2*y + 4*x^2 + 4*x*y - 4/3*x^3 - 2*x^2*y;"
674  "4/3*x^3 - 2*x^2*y;"
675  "4/3 - 4*x - 2*y + 4*x^2 + 4*x*y - 4*x^2*y + 2/3*y^3;"
676  "-2/3*y^3;"
677  "-2/3 + 2*y + 2*x^2 - 4*x*y - 4/3*x^3 + 2*x^2*y;"
678  "-2*x^2 + 4/3*x^3 + 2*x^2*y;"
679  "-2/3 + 2*y - 2*y^2 + 2/3*y^3;"
680  "-4*x*y + 2*y^2 + 4*x^2*y - 2/3*y^3;");
681 
682  bgeot::pconvex_ref cr = bgeot::parallelepiped_of_reference(2);
683  mref_convex() = cr;
684  dim() = cr->structure()->dim();
685  is_polynomialcomp() = true;
686  is_equivalent() = false;
687  is_polynomial() = false;
688  is_lagrange() = false;
689  is_standard() = false;
690  estimated_degree() = 5;
691  init_cvs_node();
692 
693  base()=std::vector<bgeot::polynomial_composite>
694  (16, bgeot::polynomial_composite(mp, false));
695  for (size_type k = 0; k < 16; ++k)
696  for (size_type ic = 0; ic < 4; ++ic)
697  base()[k].set_poly_of_subelt(ic, bgeot::read_base_poly(2, s));
698 
699  for (size_type i = 0; i < 4; ++i) {
700  base_node pt(0.0, 0.0);
701  if (i & 1) pt[0] = 1.0;
702  if (i & 2) pt[1] = 1.0;
703  add_node(lagrange_dof(2), pt);
704  add_node(derivative_dof(2, 0), pt);
705  add_node(derivative_dof(2, 1), pt);
706  }
707 
708  add_node(normal_derivative_dof(2), base_node(1.0, 0.5));
709  add_node(normal_derivative_dof(2), base_node(0.0, 0.5));
710  add_node(normal_derivative_dof(2), base_node(0.5, 1.0));
711  add_node(normal_derivative_dof(2), base_node(0.5, 0.0));
712  }
713 
714 
715  pfem quadc1p3_fem
716  (fem_param_list &params,
717  std::vector<dal::pstatic_stored_object> &dependencies) {
718  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
719  << params.size() << " should be 0.");
720  pfem p = std::make_shared<quadc1p3__>();
721  dependencies.push_back(p->ref_convex(0));
722  dependencies.push_back(p->node_tab(0));
723  return p;
724  }
725 
726  /* ******************************************************************** */
727  /* Reduced C1 composite element on quadrilateral (piecewise P3). */
728  /* ******************************************************************** */
729 
730  struct reduced_quadc1p3__ : public fem<bgeot::polynomial_composite> {
731  const quadc1p3__ *HCT;
732  virtual void mat_trans(base_matrix &M, const base_matrix &G,
733  bgeot::pgeometric_trans pgt) const;
734  virtual size_type nb_base(size_type) const { return 16; }
735  mutable base_matrix P, Mhct;
736  reduced_quadc1p3__(void);
737  };
738 
739  void reduced_quadc1p3__::mat_trans(base_matrix &M, const base_matrix &G,
740  bgeot::pgeometric_trans pgt) const {
741  HCT->mat_trans(Mhct, G, pgt);
742 
743  P(13, 1)=HCT->true_normals[1][0]*0.5; P(15, 1)=HCT->true_normals[3][0]*0.5;
744  P(13, 2)=HCT->true_normals[1][1]*0.5; P(15, 2)=HCT->true_normals[3][1]*0.5;
745 
746  P(12, 4)=HCT->true_normals[0][0]*0.5; P(15, 4)=HCT->true_normals[3][0]*0.5;
747  P(12, 5)=HCT->true_normals[0][1]*0.5; P(15, 5)=HCT->true_normals[3][1]*0.5;
748 
749  P(13, 7)=HCT->true_normals[1][0]*0.5; P(14, 7)=HCT->true_normals[2][0]*0.5;
750  P(13, 8)=HCT->true_normals[1][1]*0.5; P(14, 8)=HCT->true_normals[2][1]*0.5;
751 
752  P(12,10)=HCT->true_normals[0][0]*0.5; P(14,10)=HCT->true_normals[2][0]*0.5;
753  P(12,11)=HCT->true_normals[0][1]*0.5; P(14,11)=HCT->true_normals[2][1]*0.5;
754 
755  gmm::mult(gmm::transposed(P), Mhct, M);
756  }
757 
758  reduced_quadc1p3__::reduced_quadc1p3__(void)
759  : P(16, 12), Mhct(16, 16) {
760  HCT = dynamic_cast<const quadc1p3__ *>
761  (&(*fem_descriptor("FEM_QUADC1_COMPOSITE")));
762 
763  bgeot::pconvex_ref cr = bgeot::parallelepiped_of_reference(2);
764  mref_convex() = cr;
765  dim() = cr->structure()->dim();
766  is_polynomialcomp() = true;
767  is_equivalent() = false;
768  is_polynomial() = false;
769  is_lagrange() = false;
770  is_standard() = false;
771  estimated_degree() = 5;
772  base() = HCT->base();
773 
774  gmm::copy(gmm::identity_matrix(), P);
775  init_cvs_node();
776 
777  for (size_type i = 0; i < 4; ++i) {
778  base_node pt(0.0, 0.0);
779  if (i & 1) pt[0] = 1.0;
780  if (i & 2) pt[1] = 1.0;
781  add_node(lagrange_dof(2), pt);
782  add_node(derivative_dof(2, 0), pt);
783  add_node(derivative_dof(2, 1), pt);
784  }
785  }
786 
787 
788  pfem reduced_quadc1p3_fem
789  (fem_param_list &params,
790  std::vector<dal::pstatic_stored_object> &dependencies) {
791  GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
792  << params.size() << " should be 0.");
793  pfem p = std::make_shared<reduced_quadc1p3__>();
794  dependencies.push_back(p->ref_convex(0));
795  dependencies.push_back(p->node_tab(0));
796  return p;
797  }
798 
799 
800  /* ******************************************************************** */
801  /* HHO methods: First method for the interior of the elements and */
802  /* a method for each face (or a single method for all faces) */
803  /* ******************************************************************** */
804  /* It is guaranted (and used) that the sub-element of index 0 is the */
805  /* element itself and the faces follows in their usual order. */
806  /* It has also to be guaranted that the internal degrees of freedom are */
807  /* first. This is ensred by the dof enumeration of mesh_fem object */
808  /* since the interior element has the index 0. */
809 
810  pfem hho_method(fem_param_list &params,
811  std::vector<dal::pstatic_stored_object> &dependencies) {
812  GMM_ASSERT1(params.size() >= 2, "Bad number of parameters : "
813  << params.size() << " should be at least 2.");
814  GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
815  "Bad type of parameters");
816  pfem pf = params[0].method();
817 
818  size_type nbf = pf->ref_convex(0)->structure()->nb_faces();
819  GMM_ASSERT1(pf->is_polynomial(), "Only for polynomial elements");
820 
821  std::vector<pfem> pff(nbf);
822  if (params.size() == 2)
823  std::fill(pff.begin(), pff.end(), params[1].method());
824  else {
825  GMM_ASSERT1(params.size() == nbf+1, "Bad number of parameters : "
826  << params.size() << " a single method for all the faces or "
827  " a method for each face.");
828  for (size_type i = 0; i < nbf; ++i) {
829  GMM_ASSERT1(params[i+1].type() == 1, "Bad type of parameters");
830  GMM_ASSERT1(params[i+1].method()->is_polynomial(),
831  "Only for polynomial elements");
832  pff[i] = params[i+1].method();
833  }
834  }
835 
836  // Obtain the reference element
837  bgeot::pbasic_mesh pm;
838  bgeot::pmesh_precomposite pmp;
839  structured_mesh_for_convex(pf->ref_convex(0), 1, pm, pmp);
840 
841  // Addition of faces
842  bgeot::basic_mesh m0(*pm);
843  for (short_type i = 0; i < nbf; ++i) {
844  bgeot::mesh_structure::ind_pt_face_ct
845  ipts=m0.ind_points_of_face_of_convex(0,i);
846  bgeot::pconvex_structure cvs
847  = m0.structure_of_convex(0)->faces_structure()[i];
848  m0.add_convex(bgeot::default_trans_of_cvs(cvs), ipts.begin());
849  }
850 
851  // Build the mesh_fem
852  mesh m1(m0);
853  bgeot::mesh_precomposite mp; mp.initialise(m1);
854  mesh_fem mf(m1);
855  mf.set_finite_element(0, pf);
856  for (size_type i = 0; i < nbf; ++i)
857  mf.set_finite_element(i+1, pff[i]);
858 
859  pfem p = composite_fe_method(m1, mf, pf->ref_convex(0), true);
860  dependencies.push_back(p->ref_convex(0));
861  dependencies.push_back(p->node_tab(0));
862  return p;
863  }
864 
865  pfem interior_fem_of_hho_method(pfem hho_method) {
866 
867  const polynomial_composite_fem *phho
868  = dynamic_cast<const polynomial_composite_fem*>(hho_method.get());
869 
870  if (phho) {
871  pfem pf0 = phho->mf.fem_of_element(0);
872  pfem pf1 = phho->mf.fem_of_element(1);
873  if (pf1 && (pf1->dim()+1 == pf0->dim()))
874  return phho->mf.fem_of_element(0);
875  }
876 
877  // GMM_WARNING2("probably not a HHO method");
878  return hho_method;
879  }
880 
881 
882 
883 
884 
885 
886 
887 
888 
889 
890 
891 } /* end of namespace getfem. */
bgeot_poly_composite.h
Handle composite polynomials.
bgeot::small_vector< scalar_type >
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:98
getfem::mesh::add_triangle
size_type add_triangle(size_type a, size_type b, size_type c)
Add a triangle to the mesh, given the point id of its vertices.
Definition: getfem_mesh.cc:287
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:272
gmm::tab_ref_index_ref
indexed array reference (given a container X, and a set of indexes I, this class provides a pseudo-co...
Definition: gmm_ref.h:303
dal_naming_system.h
Naming system.
getfem_integration.h
Integration methods (exact and approximated) on convexes.
getfem_mesh_fem.h
Define the getfem::mesh_fem class.
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::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
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::derivative_dof
pdof_description derivative_dof(dim_type d, dim_type r)
Description of a unique dof of derivative type.
Definition: getfem_fem.cc:487
getfem::lagrange_dof
pdof_description lagrange_dof(dim_type d)
Description of a unique dof of lagrange type (value at the node).
Definition: getfem_fem.cc:390
getfem::normal_derivative_dof
pdof_description normal_derivative_dof(dim_type d)
Description of a unique dof of normal derivative type (normal derivative at the node,...
Definition: getfem_fem.cc:507
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
getfem::ppolycompfem
const fem< bgeot::polynomial_composite > * ppolycompfem
Polynomial composite FEM.
Definition: getfem_fem.h:600
getfem::fem_descriptor
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
Definition: getfem_fem.cc:4659
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
getfem::fem_precomp
pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt, dal::pstatic_stored_object dep)
Handles precomputations for FEM.
Definition: getfem_fem.cc:4760
getfem::ppolyfem
const fem< bgeot::base_poly > * ppolyfem
Classical polynomial FEM.
Definition: getfem_fem.h:598
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::read_base_poly
base_poly read_base_poly(short_type n, std::istream &f)
read a base_poly on the stream ist.
Definition: bgeot_poly.cc:210
bgeot::simplex_of_reference
pconvex_ref simplex_of_reference(dim_type nc, short_type K)
returns a simplex of reference of dimension nc and degree k
Definition: bgeot_convex_ref.cc:330
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:53
bgeot::structured_mesh_for_convex
void structured_mesh_for_convex(pconvex_ref cvr, short_type k, pbasic_mesh &pm, pmesh_precomposite &pmp, bool force_simplexification)
This function returns a mesh in pm which contains a refinement of the convex cvr if force_simplexific...
Definition: bgeot_poly_composite.cc:498
bgeot::parallelepiped_of_reference
pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k)
parallelepiped of reference of dimension nc (and degree 1)
Definition: bgeot_convex_ref.cc:737
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