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