GetFEM  5.4.4
getfem_fourth_order.cc
1 /*===========================================================================
2 
3  Copyright (C) 2009-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_fourth_order.h"
24 #include "getfem/getfem_omp.h"
25 
26 
27 namespace getfem {
28 
29  // ----------------------------------------------------------------------
30  // Bilaplacian brick
31  // ----------------------------------------------------------------------
32 
33  struct bilap_brick : public virtual_brick {
34 
35  virtual void asm_real_tangent_terms(const model &md, size_type ib,
36  const model::varnamelist &vl,
37  const model::varnamelist &dl,
38  const model::mimlist &mims,
39  model::real_matlist &matl,
40  model::real_veclist &,
41  model::real_veclist &,
42  size_type region,
43  build_version version) const {
44  GMM_ASSERT1(matl.size() == 1,
45  "Bilaplacian brick has one and only one term");
46  GMM_ASSERT1(mims.size() == 1,
47  "Bilaplacian brick need one and only one mesh_im");
48  GMM_ASSERT1(vl.size() == 1 && (dl.size() == 1 || dl.size() == 2),
49  "Wrong number of variables for bilaplacian brick");
50 
51  bool KL = (dl.size() == 2);
52 
53  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
54  || md.is_var_newer_than_brick(dl[0], ib);
55 
56 
57  if (recompute_matrix) {
58 
59  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
60  const mesh &m = mf_u.linked_mesh();
61  size_type Q = mf_u.get_qdim();
62  GMM_ASSERT1(Q == 1, "Bilaplacian brick is only for a scalar field");
63  const mesh_im &mim = *mims[0];
64  mesh_region rg(region);
65  m.intersect_with_mpi_region(rg);
66  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
67  const model_real_plain_vector *data = &(md.real_variable(dl[0]));
68  size_type sl = gmm::vect_size(*data);
69  if (mf_data) sl = sl * mf_data->get_qdim() / mf_data->nb_dof();
70 
71  GMM_ASSERT1(sl == 1, "Bad format of bilaplacian coefficient");
72 
73  const mesh_fem *mf_data2 = 0;
74  const model_real_plain_vector *data2 = 0;
75  if (KL) {
76  mf_data2 = md.pmesh_fem_of_variable(dl[1]);
77  data2 = &(md.real_variable(dl[1]));
78  size_type sl2 = gmm::vect_size(*data2);
79  if (mf_data2) sl = sl * mf_data2->get_qdim() / mf_data2->nb_dof();
80  GMM_ASSERT1(sl2 == 1, "Bad format of bilaplacian coefficient");
81  }
82 
83  if (KL) {
84  GMM_TRACE2("Stiffness matrix assembly of a bilaplacian term for a "
85  "Kirchhoff-Love plate");
86  }
87  else {
88  GMM_TRACE2("Stiffness matrix assembly of a bilaplacian term");
89  }
90 
91  gmm::clear(matl[0]);
92  if (mf_data) {
93  if (KL)
94  asm_stiffness_matrix_for_bilaplacian_KL
95  (matl[0], mim, mf_u, *mf_data, *data, *data2, rg);
96  else
97  asm_stiffness_matrix_for_bilaplacian
98  (matl[0], mim, mf_u, *mf_data, *data, rg);
99  } else {
100  if (KL) {
101  asm_stiffness_matrix_for_homogeneous_bilaplacian_KL
102  (matl[0], mim, mf_u, *data, *data2, rg);
103  }
104  else
105  asm_stiffness_matrix_for_homogeneous_bilaplacian
106  (matl[0], mim, mf_u, *data, rg);
107  }
108  }
109  }
110 
111  bilap_brick(void) {
112  set_flags("Bilaplacian operator", true /* is linear*/,
113  true /* is symmetric */, true /* is coercive */,
114  true /* is real */, false /* is complex */);
115  }
116 
117  };
118 
119  size_type add_bilaplacian_brick
120  (model &md, const mesh_im &mim, const std::string &varname,
121  const std::string &dataname,
122  size_type region) {
123  pbrick pbr = std::make_shared<bilap_brick>();
124  model::termlist tl;
125  tl.push_back(model::term_description(varname, varname, true));
126  model::varnamelist dl(1, dataname);
127  return md.add_brick(pbr, model::varnamelist(1, varname), dl, tl,
128  model::mimlist(1, &mim), region);
129  }
130 
131  size_type add_bilaplacian_brick_KL
132  (model &md, const mesh_im &mim, const std::string &varname,
133  const std::string &dataname1, const std::string &dataname2,
134  size_type region) {
135  pbrick pbr = std::make_shared<bilap_brick>();
136  model::termlist tl;
137  tl.push_back(model::term_description(varname, varname, true));
138  model::varnamelist dl(1, dataname1);
139  dl.push_back(dataname2);
140  return md.add_brick(pbr, model::varnamelist(1, varname), dl, tl,
141  model::mimlist(1, &mim), region);
142  }
143 
144 
145 
146  // ----------------------------------------------------------------------
147  //
148  // Normal derivative source term brick
149  //
150  // ----------------------------------------------------------------------
151 
152  struct normal_derivative_source_term_brick : public virtual_brick {
153 
154  virtual void asm_real_tangent_terms(const model &md, size_type,
155  const model::varnamelist &vl,
156  const model::varnamelist &dl,
157  const model::mimlist &mims,
158  model::real_matlist &,
159  model::real_veclist &vecl,
160  model::real_veclist &,
161  size_type region,
162  build_version) const {
163  GMM_ASSERT1(vecl.size() == 1,
164  "Normal derivative source term brick has one and only "
165  "one term");
166  GMM_ASSERT1(mims.size() == 1,
167  "Normal derivative source term brick need one and only "
168  "one mesh_im");
169  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
170  "Wrong number of variables for normal derivative "
171  "source term brick");
172 
173  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
174  const mesh_im &mim = *mims[0];
175  const model_real_plain_vector &A = md.real_variable(dl[0]);
176  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
177  mesh_region rg(region);
178  mim.linked_mesh().intersect_with_mpi_region(rg);
179 
180  size_type s = gmm::vect_size(A);
181  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
182 
183  GMM_ASSERT1(s == mf_u.get_qdim()
184  || s == size_type(mf_u.get_qdim()*gmm::sqr(mf_u.linked_mesh().dim())),
185  dl[0] << ": bad format of normal derivative source term "
186  "data. Detected dimension is " << s << " should be "
187  << size_type(mf_u.get_qdim()));
188 
189  GMM_TRACE2("Normal derivative source term assembly");
190  if (mf_data)
191  asm_normal_derivative_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
192  else
193  asm_homogeneous_normal_derivative_source_term(vecl[0], mim, mf_u, A, rg);
194 
195  }
196 
197  virtual void asm_complex_tangent_terms(const model &md, size_type,
198  const model::varnamelist &vl,
199  const model::varnamelist &dl,
200  const model::mimlist &mims,
201  model::complex_matlist &,
202  model::complex_veclist &vecl,
203  model::complex_veclist &,
204  size_type region,
205  build_version) const {
206  GMM_ASSERT1(vecl.size() == 1,
207  "Normal derivative source term brick has one and only "
208  "one term");
209  GMM_ASSERT1(mims.size() == 1,
210  "Normal derivative source term brick need one and only "
211  "one mesh_im");
212  GMM_ASSERT1(vl.size() == 1 && dl.size() == 1,
213  "Wrong number of variables for normal derivative "
214  "source term brick");
215 
216  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
217  const mesh_im &mim = *mims[0];
218  const model_complex_plain_vector &A = md.complex_variable(dl[0]);
219  const mesh_fem *mf_data = md.pmesh_fem_of_variable(dl[0]);
220  mesh_region rg(region);
221  mim.linked_mesh().intersect_with_mpi_region(rg);
222 
223  size_type s = gmm::vect_size(A);
224  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
225 
226  GMM_ASSERT1(mf_u.get_qdim() == s,
227  dl[0] << ": bad format of normal derivative source term "
228  "data. Detected dimension is " << s << " should be "
229  << size_type(mf_u.get_qdim()));
230 
231  GMM_TRACE2("Normal derivative source term assembly");
232  if (mf_data)
233  asm_normal_derivative_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
234  else
235  asm_homogeneous_normal_derivative_source_term(vecl[0], mim, mf_u, A, rg);
236 
237  }
238 
239  normal_derivative_source_term_brick(void) {
240  set_flags("Normal derivative source term", true /* is linear*/,
241  true /* is symmetric */, true /* is coercive */,
242  true /* is real */, true /* is complex */,
243  false /* compute each time */);
244  }
245 
246 
247  };
248 
249  size_type add_normal_derivative_source_term_brick
250  (model &md, const mesh_im &mim, const std::string &varname,
251  const std::string &dataname, size_type region) {
252  pbrick pbr = std::make_shared<normal_derivative_source_term_brick>();
253  model::termlist tl;
254  tl.push_back(model::term_description(varname));
255  model::varnamelist vdata(1, dataname);
256  return md.add_brick(pbr, model::varnamelist(1, varname),
257  vdata, tl, model::mimlist(1, &mim), region);
258  }
259 
260 
261 
262  // ----------------------------------------------------------------------
263  //
264  // K.L. source term brick
265  //
266  // ----------------------------------------------------------------------
267 
268  struct KL_source_term_brick : public virtual_brick {
269 
270  virtual void asm_real_tangent_terms(const model &md, size_type,
271  const model::varnamelist &vl,
272  const model::varnamelist &dl,
273  const model::mimlist &mims,
274  model::real_matlist &,
275  model::real_veclist &vecl,
276  model::real_veclist &,
277  size_type region,
278  build_version) const {
279  GMM_ASSERT1(vecl.size() == 1,
280  "Kirchhoff Love source term brick has one and only "
281  "one term");
282  GMM_ASSERT1(mims.size() == 1,
283  "Kirchhoff Love source term brick need one and only "
284  "one mesh_im");
285  GMM_ASSERT1(vl.size() == 1 && dl.size() == 2,
286  "Wrong number of variables for Kirchhoff Love "
287  "source term brick");
288 
289  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
290  const mesh_im &mim = *mims[0];
291  const model_real_plain_vector &A = md.real_variable(dl[0]);
292  const mesh_fem *mf_dataA = md.pmesh_fem_of_variable(dl[0]);
293  const model_real_plain_vector &B = md.real_variable(dl[1]);
294  const mesh_fem *mf_dataB = md.pmesh_fem_of_variable(dl[1]);
295  size_type N = mf_u.linked_mesh().dim();
296  mesh_region rg(region);
297  mim.linked_mesh().intersect_with_mpi_region(rg);
298 
299  size_type s = gmm::vect_size(A);
300  if (mf_dataA) s = s * mf_dataA->get_qdim() / mf_dataA->nb_dof();
301  GMM_ASSERT1(mf_u.get_qdim() == 1 && s == N*N,
302  dl[0] << ": bad format of Kirchhoff Love Neumann term "
303  "data. Detected dimension is " << s << " should be "
304  << N*N);
305 
306  s = gmm::vect_size(B);
307  if (mf_dataB) s = s * mf_dataB->get_qdim() / mf_dataB->nb_dof();
308  GMM_ASSERT1(s == N,
309  dl[0] << ": bad format of Kirchhoff Love Neumann term "
310  "data. Detected dimension is " << s << " should be "
311  << N);
312 
313 
314  GMM_TRACE2("Kirchhoff Love Neumann term assembly");
315  if (mf_dataA)
316  asm_neumann_KL_term(vecl[0], mim, mf_u, *mf_dataA, A, B, rg);
317  else
318  asm_neumann_KL_homogeneous_term(vecl[0], mim, mf_u, A, B, rg);
319 
320  }
321 
322  KL_source_term_brick(void) {
323  set_flags("Kirchhoff Love Neumann term", true /* is linear*/,
324  true /* is symmetric */, true /* is coercive */,
325  true /* is real */, false /* is complex */,
326  false /* compute each time */);
327  }
328 
329 
330  };
331 
332  size_type add_Kirchhoff_Love_Neumann_term_brick
333  (model &md, const mesh_im &mim, const std::string &varname,
334  const std::string &dataname1, const std::string &dataname2,
335  size_type region) {
336  pbrick pbr = std::make_shared<KL_source_term_brick>();
337  model::termlist tl;
338  tl.push_back(model::term_description(varname));
339  model::varnamelist vdata(1, dataname1);
340  vdata.push_back(dataname2);
341  return md.add_brick(pbr, model::varnamelist(1, varname),
342  vdata, tl, model::mimlist(1, &mim), region);
343  }
344 
345 
346 
347 
348 
349 
350 
351 
352 
353 
354 
355 
356 
357 
358 
359 
360 
361 
362 
363 
364 
365 
366 
367 
368  // ----------------------------------------------------------------------
369  //
370  // Normal derivative Dirichlet condition brick
371  //
372  // ----------------------------------------------------------------------
373  // Two variables : with multipliers
374  // One variable : penalization
375 
376  struct normal_derivative_Dirichlet_condition_brick : public virtual_brick {
377 
378  bool R_must_be_derivated;
379  mutable getfem::omp_distribute<model_real_sparse_matrix> rB_th;
380  mutable getfem::omp_distribute<model_real_plain_vector> rV_th;
381  mutable getfem::omp_distribute<model_complex_sparse_matrix> cB_th;
382  mutable getfem::omp_distribute<model_complex_plain_vector> cV_th;
383 
384  virtual void asm_real_tangent_terms(const model &md, size_type ib,
385  const model::varnamelist &vl,
386  const model::varnamelist &dl,
387  const model::mimlist &mims,
388  model::real_matlist &matl,
389  model::real_veclist &vecl,
390  model::real_veclist &,
391  size_type region,
392  build_version version) const {
393  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
394  "Normal derivative Dirichlet condition brick has one and "
395  "only one term");
396  GMM_ASSERT1(mims.size() == 1,
397  "Normal derivative Dirichlet condition brick need one and "
398  "only one mesh_im");
399  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 2,
400  "Wrong number of variables for normal derivative Dirichlet "
401  "condition brick");
402 
403  model_real_sparse_matrix& rB = rB_th;
404  model_real_plain_vector& rV = rV_th;
405 
406  bool penalized = (vl.size() == 1);
407  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
408  const mesh_fem &mf_mult = md.mesh_fem_of_variable(vl[vl.size()-1]);
409  const mesh_im &mim = *mims[0];
410  const model_real_plain_vector *A = 0, *COEFF = 0;
411  const mesh_fem *mf_data = 0;
412  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
413  || (penalized && md.is_var_newer_than_brick(dl[0], ib));
414 
415  if (penalized) {
416  COEFF = &(md.real_variable(dl[0]));
417  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
418  "Data for coefficient should be a scalar");
419  }
420 
421  size_type s = 0, ind = (penalized ? 1 : 0);
422  if (dl.size() > ind) {
423  A = &(md.real_variable(dl[ind]));
424  mf_data = md.pmesh_fem_of_variable(dl[ind]);
425  s = gmm::vect_size(*A);
426  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
427  GMM_ASSERT1(s == mf_u.get_qdim() || s == mf_u.linked_mesh().dim(),
428  dl[ind] << ": bad format of normal derivative Dirichlet "
429  "data. Detected dimension is " << s << " should be "
430  << size_type(mf_u.get_qdim()) << " or "
431  << mf_u.linked_mesh().dim());
432  }
433 
434  mesh_region rg(region);
435  mim.linked_mesh().intersect_with_mpi_region(rg);
436 
437  if (recompute_matrix) {
438  GMM_TRACE2("Mass term assembly for normal derivative Dirichlet "
439  "condition");
440  if (penalized) {
441  gmm::resize(rB, mf_mult.nb_dof(), mf_u.nb_dof());
442  gmm::clear(rB);
443  asm_normal_derivative_dirichlet_constraints
444  (rB, vecl[0], mim, mf_u, mf_mult,
445  *mf_data, *A, rg, R_must_be_derivated, ASMDIR_BUILDH);
446  } else {
447  gmm::clear(matl[0]);
448  asm_normal_derivative_dirichlet_constraints
449  (matl[0], vecl[0], mim, mf_u, mf_mult,
450  *mf_data, *A, rg, R_must_be_derivated, ASMDIR_BUILDH);
451  }
452 
453  if (penalized) {
454  gmm::mult(gmm::transposed(rB), rB, matl[0]);
455  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
456  }
457  }
458 
459  if (dl.size() > ind) {
460  GMM_TRACE2("Source term assembly for normal derivative Dirichlet "
461  "condition");
462  model_real_plain_vector *R = penalized ? &rV : &(vecl[0]);
463  if (penalized) { gmm::resize(rV, mf_mult.nb_dof()); gmm::clear(rV); }
464 
465  if (mf_data) {
466  if (!R_must_be_derivated) {
467  if (s == mf_u.linked_mesh().dim())
468  asm_normal_source_term(*R, mim, mf_mult, *mf_data, *A, rg);
469  else
470  asm_source_term(*R, mim, mf_mult, *mf_data, *A, rg);
471  }
472  else {
473  asm_real_or_complex_1_param_vec
474  (*R, mim, mf_mult, mf_data, *A, rg, "(Grad_A.Normal)*Test_u");
475  // asm_real_or_complex_1_param
476  // (*R, mim, mf_mult, *mf_data, *A, rg,
477  // "R=data(#2); V(#1)+=comp(Base(#1).Grad(#2).Normal())(:,i,j,j).R(i)");
478  }
479  } else {
480  GMM_ASSERT1(!R_must_be_derivated, "Incoherent situation");
481  if (s == mf_u.linked_mesh().dim())
482  asm_homogeneous_normal_source_term(*R, mim, mf_mult, *A, rg);
483  else
484  asm_homogeneous_source_term(*R, mim, mf_mult, *A, rg);
485  }
486  if (penalized) {
487  gmm::mult(gmm::transposed(rB), rV, vecl[0]);
488  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
489  rV = model_real_plain_vector();
490  }
491  }
492 
493 
494  }
495 
496  virtual void asm_complex_tangent_terms(const model &md, size_type ib,
497  const model::varnamelist &vl,
498  const model::varnamelist &dl,
499  const model::mimlist &mims,
500  model::complex_matlist &matl,
501  model::complex_veclist &vecl,
502  model::complex_veclist &,
503  size_type region,
504  build_version version) const {
505  GMM_ASSERT1(vecl.size() == 1 && matl.size() == 1,
506  "Normal derivative Dirichlet condition brick has one and "
507  "only one term");
508  GMM_ASSERT1(mims.size() == 1,
509  "Normal derivative Dirichlet condition brick need one and "
510  "only one mesh_im");
511  GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() <= 2,
512  "Wrong number of variables for normal derivative Dirichlet "
513  "condition brick");
514 
515  model_complex_sparse_matrix& cB = cB_th;
516  model_complex_plain_vector& cV = cV_th;
517 
518  bool penalized = (vl.size() == 1);
519  const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
520  const mesh_fem &mf_mult = md.mesh_fem_of_variable(vl[vl.size()-1]);
521  const mesh_im &mim = *mims[0];
522  const model_complex_plain_vector *A = 0, *COEFF = 0;
523  const mesh_fem *mf_data = 0;
524  bool recompute_matrix = !((version & model::BUILD_ON_DATA_CHANGE) != 0)
525  || (penalized && md.is_var_newer_than_brick(dl[0], ib));
526 
527  if (penalized) {
528  COEFF = &(md.complex_variable(dl[0]));
529  GMM_ASSERT1(gmm::vect_size(*COEFF) == 1,
530  "Data for coefficient should be a scalar");
531  }
532 
533  size_type s = 0, ind = (penalized ? 1 : 0);
534  if (dl.size() > ind) {
535  A = &(md.complex_variable(dl[ind]));
536  mf_data = md.pmesh_fem_of_variable(dl[ind]);
537  s = gmm::vect_size(*A);
538  if (mf_data) s = s * mf_data->get_qdim() / mf_data->nb_dof();
539  GMM_ASSERT1(s == mf_u.get_qdim() || s == mf_u.linked_mesh().dim(),
540  dl[ind] << ": bad format of normal derivative Dirichlet "
541  "data. Detected dimension is " << s << " should be "
542  << size_type(mf_u.get_qdim()) << " or "
543  << mf_u.linked_mesh().dim());
544  }
545 
546  mesh_region rg(region);
547  mim.linked_mesh().intersect_with_mpi_region(rg);
548 
549  if (recompute_matrix) {
550  GMM_TRACE2("Mass term assembly for normal derivative Dirichlet "
551  "condition");
552  if (penalized) {
553  gmm::resize(cB, mf_mult.nb_dof(), mf_u.nb_dof());
554  gmm::clear(cB);
555  asm_normal_derivative_dirichlet_constraints
556  (cB, vecl[0], mim, mf_u, mf_mult,
557  *mf_data, *A, rg, R_must_be_derivated, ASMDIR_BUILDH);
558  } else {
559  gmm::clear(matl[0]);
560  asm_normal_derivative_dirichlet_constraints
561  (matl[0], vecl[0], mim, mf_u, mf_mult,
562  *mf_data, *A, rg, R_must_be_derivated, ASMDIR_BUILDH);
563  }
564 
565  if (penalized) {
566  gmm::mult(gmm::transposed(cB), cB, matl[0]);
567  gmm::scale(matl[0], gmm::abs((*COEFF)[0]));
568  }
569  }
570 
571  if (dl.size() > ind) {
572  GMM_TRACE2("Source term assembly for normal derivative Dirichlet "
573  "condition");
574  model_complex_plain_vector *R = penalized ? &cV : &(vecl[0]);
575  if (penalized) { gmm::resize(cV, mf_mult.nb_dof()); gmm::clear(cV); }
576 
577  if (mf_data) {
578  if (!R_must_be_derivated) {
579  if (s == mf_u.linked_mesh().dim())
580  asm_normal_source_term(*R, mim, mf_mult, *mf_data, *A, rg);
581  else
582  asm_source_term(*R, mim, mf_mult, *mf_data, *A, rg);
583  }
584  else {
585  asm_real_or_complex_1_param_vec
586  (*R, mim, mf_mult, mf_data, *A, rg, "(Grad_A.Normal)*Test_u");
587  // asm_real_or_complex_1_param
588  // (*R, mim, mf_mult, *mf_data, *A, rg,
589  // "R=data(#2); V(#1)+=comp(Base(#1).Grad(#2).Normal())(:,i,j,j).R(i)");
590  }
591  } else {
592  GMM_ASSERT1(!R_must_be_derivated, "Incoherent situation");
593  if (s == mf_u.linked_mesh().dim())
594  asm_homogeneous_normal_source_term(*R, mim, mf_mult, *A, rg);
595  else
596  asm_homogeneous_source_term(*R, mim, mf_mult, *A, rg);
597  }
598  if (penalized) {
599  gmm::mult(gmm::transposed(cB), cV, vecl[0]);
600  gmm::scale(vecl[0], gmm::abs((*COEFF)[0]));
601  cV = model_complex_plain_vector();
602  }
603  }
604  }
605 
606  normal_derivative_Dirichlet_condition_brick
607  (bool penalized, bool R_must_be_derivated_) {
608  R_must_be_derivated = R_must_be_derivated_;
609  set_flags(penalized ?
610  "Normal derivative Dirichlet with penalization brick"
611  : "Normal derivative Dirichlet with multipliers brick",
612  true /* is linear*/,
613  true /* is symmetric */, penalized /* is coercive */,
614  true /* is real */, true /* is complex */,
615  false /* compute each time */);
616  }
617  };
618 
619  size_type add_normal_derivative_Dirichlet_condition_with_multipliers
620  (model &md, const mesh_im &mim, const std::string &varname,
621  const std::string &multname, size_type region,
622  const std::string &dataname, bool R_must_be_derivated) {
623  pbrick pbr = std::make_shared<normal_derivative_Dirichlet_condition_brick>(false, R_must_be_derivated);
624  model::termlist tl;
625  tl.push_back(model::term_description(multname, varname, true));
626  model::varnamelist vl(1, varname);
627  vl.push_back(multname);
628  model::varnamelist dl;
629  if (dataname.size()) dl.push_back(dataname);
630  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
631  }
632 
633  size_type add_normal_derivative_Dirichlet_condition_with_multipliers
634  (model &md, const mesh_im &mim, const std::string &varname,
635  const mesh_fem &mf_mult, size_type region,
636  const std::string &dataname, bool R_must_be_derivated) {
637  std::string multname = md.new_name("mult_on_" + varname);
638  md.add_multiplier(multname, mf_mult, varname);
639  return add_normal_derivative_Dirichlet_condition_with_multipliers
640  (md, mim, varname, multname, region, dataname, R_must_be_derivated);
641  }
642 
643  size_type add_normal_derivative_Dirichlet_condition_with_multipliers
644  (model &md, const mesh_im &mim, const std::string &varname,
645  dim_type degree, size_type region,
646  const std::string &dataname, bool R_must_be_derivated) {
647  const mesh_fem &mf_u = md.mesh_fem_of_variable(varname);
648  const mesh_fem &mf_mult = classical_mesh_fem(mf_u.linked_mesh(),
649  degree, mf_u.get_qdim());
650  return add_normal_derivative_Dirichlet_condition_with_multipliers
651  (md, mim, varname, mf_mult, region, dataname, R_must_be_derivated);
652  }
653 
654 
655  size_type add_normal_derivative_Dirichlet_condition_with_penalization
656  (model &md, const mesh_im &mim, const std::string &varname,
657  scalar_type penalisation_coeff, size_type region,
658  const std::string &dataname, bool R_must_be_derivated) {
659  std::string coeffname = md.new_name("penalization_on_" + varname);
660  md.add_fixed_size_data(coeffname, 1);
661  if (md.is_complex())
662  md.set_complex_variable(coeffname)[0] = penalisation_coeff;
663  else
664  md.set_real_variable(coeffname)[0] = penalisation_coeff;
665  pbrick pbr = std::make_shared<normal_derivative_Dirichlet_condition_brick>(true, R_must_be_derivated);
666  model::termlist tl;
667  tl.push_back(model::term_description(varname, varname, true));
668  model::varnamelist vl(1, varname);
669  model::varnamelist dl(1, coeffname);
670  if (dataname.size()) dl.push_back(dataname);
671  return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
672  }
673 
674  // + changement du coeff de penalisation avec la fonction de Dirichelt standard.
675 
676 
677 
678 } /* end of namespace getfem. */
679 
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:152
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:316
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:283
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::model
`‘Model’' variables store the variables, the data and the description of a model.
Definition: getfem_models.h:115
getfem::model::add_fixed_size_data
void add_fixed_size_data(const std::string &name, size_type size, size_type niter=1)
Add a fixed size data to the model.
Definition: getfem_models.cc:762
getfem::model::add_brick
size_type add_brick(pbrick pbr, const varnamelist &varnames, const varnamelist &datanames, const termlist &terms, const mimlist &mims, size_type region)
Add a brick to the model.
Definition: getfem_models.cc:1032
getfem::model::mesh_fem_of_variable
const mesh_fem & mesh_fem_of_variable(const std::string &name) const
Gives the access to the mesh_fem of a variable if any.
Definition: getfem_models.cc:2889
getfem::model::add_multiplier
void add_multiplier(const std::string &name, const mesh_fem &mf, const std::string &primal_name, size_type niter=1)
Add a particular variable linked to a fem being a multiplier with respect to a primal variable.
Definition: getfem_models.cc:888
getfem::model::set_real_variable
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
Definition: getfem_models.cc:3011
getfem::model::is_complex
bool is_complex() const
Boolean which says if the model deals with real or complex unknowns and data.
Definition: getfem_models.h:563
getfem::model::set_complex_variable
model_complex_plain_vector & set_complex_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
Definition: getfem_models.cc:3039
getfem::model::new_name
std::string new_name(const std::string &name)
Gives a non already existing variable name begining by name.
Definition: getfem_models.cc:175
getfem::omp_distribute< model_real_sparse_matrix >
getfem_fourth_order.h
assembly procedures and bricks for fourth order pdes.
getfem_omp.h
Tools for multithreaded, OpenMP and Boost based parallelization.
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
getfem::asm_stiffness_matrix_for_bilaplacian
void asm_stiffness_matrix_for_bilaplacian(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
assembly of .
Definition: getfem_fourth_order.h:56
getfem::asm_source_term
void asm_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg=mesh_region::all_convexes())
source term (for both volumic sources and boundary (Neumann) sources).
Definition: getfem_assembling.h:877
getfem::asm_homogeneous_normal_source_term
void asm_homogeneous_normal_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const VECT2 &F, const mesh_region &rg)
Homogeneous normal source term (for boundary (Neumann) condition).
Definition: getfem_assembling.h:919
getfem::asm_normal_derivative_dirichlet_constraints
void asm_normal_derivative_dirichlet_constraints(MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_mult, const mesh_fem &mf_r, const VECT2 &r_data, const mesh_region &rg, bool R_must_be_derivated, int version)
Assembly of normal derivative Dirichlet constraints in a weak form.
Definition: getfem_fourth_order.h:352
getfem::asm_normal_derivative_source_term
void asm_normal_derivative_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg)
assembly of .
Definition: getfem_fourth_order.h:161
getfem::asm_normal_source_term
void asm_normal_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg)
Normal source term (for boundary (Neumann) condition).
Definition: getfem_assembling.h:905
getfem::asm_homogeneous_source_term
void asm_homogeneous_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const VECT2 &F, const mesh_region &rg=mesh_region::all_convexes())
source term (for both volumic sources and boundary (Neumann) sources).
Definition: getfem_assembling.h:893
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:47
getfem::add_bilaplacian_brick
size_type add_bilaplacian_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname, size_type region=size_type(-1))
Adds a bilaplacian brick on the variable varname and on the mesh region region.
Definition: getfem_fourth_order.cc:120
getfem::add_bilaplacian_brick_KL
size_type add_bilaplacian_brick_KL(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname1, const std::string &dataname2, size_type region=size_type(-1))
Adds a bilaplacian brick on the variable varname and on the mesh region region.
Definition: getfem_fourth_order.cc:132
getfem::add_normal_derivative_Dirichlet_condition_with_multipliers
size_type add_normal_derivative_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname=std::string(), bool R_must_be_derivated=false)
Adds a Dirichlet condition on the normal derivative of the variable varname and on the mesh region re...
Definition: getfem_fourth_order.cc:620
getfem::pbrick
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
Definition: getfem_models.h:49
getfem::classical_mesh_fem
const mesh_fem & classical_mesh_fem(const mesh &mesh, dim_type degree, dim_type qdim=1, bool complete=false)
Gives the descriptor of a classical finite element method of degree K on mesh.
Definition: getfem_mesh_fem.cc:861
getfem::add_normal_derivative_source_term_brick
size_type add_normal_derivative_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname, size_type region)
Adds a normal derivative source term brick :math:F = \int b.
Definition: getfem_fourth_order.cc:250
getfem::add_Kirchhoff_Love_Neumann_term_brick
size_type add_Kirchhoff_Love_Neumann_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname1, const std::string &dataname2, size_type region)
Adds a Neumann term brick for Kirchhoff-Love model on the variable varname and the mesh region region...
Definition: getfem_fourth_order.cc:333
getfem::add_normal_derivative_Dirichlet_condition_with_penalization
size_type add_normal_derivative_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalisation_coeff, size_type region, const std::string &dataname=std::string(), bool R_must_be_derivated=false)
Adds a Dirichlet condition on the normal derivative of the variable varname and on the mesh region re...
Definition: getfem_fourth_order.cc:656