GetFEM++  5.3
getfem_assembling.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-2017 Yves Renard, Julien Pommier
5 
6  This file is a part of GetFEM++
7 
8  GetFEM++ is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /** @file getfem_assembling.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>
34  @author Julien Pommier <Julien.Pommier@insa-toulouse.fr>
35  @date November 17, 2000.
36  @brief Miscelleanous assembly routines for common terms. Use the low-level
37  generic assembly. Prefer the high-level one.
38  */
39 
40 /** @defgroup asm Assembly routines */
41 
42 #ifndef GETFEM_ASSEMBLING_H__
43 #define GETFEM_ASSEMBLING_H__
44 
45 #include "getfem_assembling_tensors.h"
46 #include "getfem_generic_assembly.h"
47 
48 namespace getfem {
49 
50  /**
51  compute @f$ \|U\|_2 @f$, U might be real or complex
52  @ingroup asm
53  */
54  template<typename VEC>
55  inline scalar_type asm_L2_norm
56  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
57  const mesh_region &rg=mesh_region::all_convexes())
58  { return sqrt(asm_L2_norm_sqr(mim, mf, U, rg)); }
59 
60  template<typename VEC>
61  scalar_type asm_L2_norm_sqr
62  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
63  const mesh_region &rg=mesh_region::all_convexes()) {
64  return asm_L2_norm_sqr(mim, mf, U, rg,
65  typename gmm::linalg_traits<VEC>::value_type());
66  }
67 
68  template<typename VEC, typename T>
69  inline scalar_type asm_L2_norm_sqr(const mesh_im &mim, const mesh_fem &mf,
70  const VEC &U, const mesh_region &rg_, T) {
71  ga_workspace workspace;
72  model_real_plain_vector UU(mf.nb_dof());
73  gmm::copy(U, UU);
74  gmm::sub_interval Iu(0, mf.nb_dof());
75  workspace.add_fem_variable("u", mf, Iu, UU);
76  workspace.add_expression("u.u", mim, rg_);
77  workspace.assembly(0);
78  return workspace.assembled_potential();
79  }
80 
81  inline scalar_type asm_L2_norm_sqr(const mesh_im &mim, const mesh_fem &mf,
82  const model_real_plain_vector &U,
83  const mesh_region &rg_, scalar_type) {
84  ga_workspace workspace;
85  gmm::sub_interval Iu(0, mf.nb_dof());
86  workspace.add_fem_variable("u", mf, Iu, U);
87  workspace.add_expression("u.u", mim, rg_);
88  workspace.assembly(0);
89  return workspace.assembled_potential();
90  }
91 
92  template<typename VEC, typename T>
93  inline scalar_type asm_L2_norm_sqr(const mesh_im &mim, const mesh_fem &mf,
94  const VEC &U,
95  const mesh_region &rg, std::complex<T>) {
96  ga_workspace workspace;
97  model_real_plain_vector UUR(mf.nb_dof()), UUI(mf.nb_dof());
98  gmm::copy(gmm::real_part(U), UUR);
99  gmm::copy(gmm::imag_part(U), UUI);
100  gmm::sub_interval Iur(0, mf.nb_dof()), Iui(mf.nb_dof(), mf.nb_dof());
101  workspace.add_fem_variable("u", mf, Iur, UUR);
102  workspace.add_fem_variable("v", mf, Iui, UUI);
103  workspace.add_expression("u.u + v.v", mim, rg);
104  workspace.assembly(0);
105  return workspace.assembled_potential();
106  }
107 
108  /**
109  Compute the distance between U1 and U2, defined on two different
110  mesh_fems (but sharing the same mesh), without interpolating U1 on mf2.
111 
112  @ingroup asm
113  */
114  template<typename VEC1, typename VEC2>
115  inline scalar_type asm_L2_dist
116  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
117  const mesh_fem &mf2, const VEC2 &U2,
118  mesh_region rg = mesh_region::all_convexes()) {
119  return sqrt(asm_L2_dist_sqr(mim, mf1, U1, mf2, U2, rg,
120  typename gmm::linalg_traits<VEC1>::value_type()));
121  }
122 
123  template<typename VEC1, typename VEC2, typename T>
124  inline scalar_type asm_L2_dist_sqr
125  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
126  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, T) {
127  ga_workspace workspace;
128  model_real_plain_vector UU1(mf1.nb_dof()), UU2(mf2.nb_dof());
129  gmm::copy(U1, UU1); gmm::copy(U2, UU2);
130  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
131  workspace.add_fem_variable("u1", mf1, Iu1, UU1);
132  workspace.add_fem_variable("u2", mf2, Iu2, UU2);
133  workspace.add_expression("(u2-u1).(u2-u1)", mim, rg);
134  workspace.assembly(0);
135  return workspace.assembled_potential();
136  }
137 
138  inline scalar_type asm_L2_dist_sqr
139  (const mesh_im &mim, const mesh_fem &mf1, const model_real_plain_vector &U1,
140  const mesh_fem &mf2, const model_real_plain_vector &U2,
141  mesh_region rg, scalar_type) {
142  ga_workspace workspace;
143  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
144  workspace.add_fem_variable("u1", mf1, Iu1, U1);
145  workspace.add_fem_variable("u2", mf2, Iu2, U2);
146  workspace.add_expression("(u2-u1).(u2-u1)", mim, rg);
147  workspace.assembly(0);
148  return workspace.assembled_potential();
149  }
150 
151  template<typename VEC1, typename VEC2, typename T>
152  inline scalar_type asm_L2_dist_sqr
153  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
154  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, std::complex<T>) {
155  ga_workspace workspace;
156  model_real_plain_vector UU1R(mf1.nb_dof()), UU2R(mf2.nb_dof());
157  model_real_plain_vector UU1I(mf1.nb_dof()), UU2I(mf2.nb_dof());
158  gmm::copy(gmm::real_part(U1), UU1R); gmm::copy(gmm::imag_part(U1), UU1I);
159  gmm::copy(gmm::real_part(U2), UU2R); gmm::copy(gmm::imag_part(U2), UU2I);
160  gmm::sub_interval Iu1r(0, mf1.nb_dof()), Iu2r(mf1.nb_dof(), mf2.nb_dof());
161  gmm::sub_interval Iu1i(Iu2r.last(), mf1.nb_dof());
162  gmm::sub_interval Iu2i(Iu1i.last(), mf2.nb_dof());
163  workspace.add_fem_variable("u1", mf1, Iu1r, UU1R);
164  workspace.add_fem_variable("u2", mf2, Iu2r, UU2R);
165  workspace.add_fem_variable("v1", mf1, Iu1i, UU1I);
166  workspace.add_fem_variable("v2", mf2, Iu2i, UU2I);
167  workspace.add_expression("(u2-u1).(u2-u1) + (v2-v1).(v2-v1)", mim, rg);
168  workspace.assembly(0);
169  return workspace.assembled_potential();
170  }
171 
172 
173  /**
174  compute @f$\|\nabla U\|_2@f$, U might be real or complex
175  @ingroup asm
176  */
177  template<typename VEC>
178  scalar_type asm_H1_semi_norm
179  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
180  const mesh_region &rg = mesh_region::all_convexes()) {
181  typedef typename gmm::linalg_traits<VEC>::value_type T;
182  return sqrt(asm_H1_semi_norm_sqr(mim, mf, U, rg, T()));
183  }
184 
185  template<typename VEC, typename T>
186  inline scalar_type asm_H1_semi_norm_sqr
187  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
188  const mesh_region &rg_, T) {
189  ga_workspace workspace;
190  model_real_plain_vector UU(mf.nb_dof()); gmm::copy(U, UU);
191  gmm::sub_interval Iu(0, mf.nb_dof());
192  workspace.add_fem_variable("u", mf, Iu, UU);
193  workspace.add_expression("Grad_u:Grad_u", mim, rg_);
194  workspace.assembly(0);
195  return workspace.assembled_potential();
196  }
197 
198  inline scalar_type asm_H1_semi_norm_sqr
199  (const mesh_im &mim, const mesh_fem &mf, const model_real_plain_vector &U,
200  const mesh_region &rg_, scalar_type) {
201  ga_workspace workspace;
202  gmm::sub_interval Iu(0, mf.nb_dof());
203  workspace.add_fem_variable("u", mf, Iu, U);
204  workspace.add_expression("Grad_u:Grad_u", mim, rg_);
205  workspace.assembly(0);
206  return workspace.assembled_potential();
207  }
208 
209 
210  template<typename VEC, typename T>
211  inline scalar_type asm_H1_semi_norm_sqr
212  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
213  const mesh_region &rg, std::complex<T>) {
214  ga_workspace workspace;
215  model_real_plain_vector UUR(mf.nb_dof()), UUI(mf.nb_dof());
216  gmm::copy(gmm::real_part(U), UUR);
217  gmm::copy(gmm::imag_part(U), UUI);
218  gmm::sub_interval Iur(0, mf.nb_dof()), Iui(mf.nb_dof(), mf.nb_dof());
219  workspace.add_fem_variable("u", mf, Iur, UUR);
220  workspace.add_fem_variable("v", mf, Iui, UUI);
221  workspace.add_expression("Grad_u:Grad_u + Grad_v:Grad_v", mim, rg);
222  workspace.assembly(0);
223  return workspace.assembled_potential();
224  }
225 
226  /**
227  Compute the H1 semi-distance between U1 and U2, defined on two different
228  mesh_fems (but sharing the same mesh), without interpolating U1 on mf2.
229 
230  @ingroup asm
231  */
232  template<typename VEC1, typename VEC2>
233  inline scalar_type asm_H1_semi_dist
234  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
235  const mesh_fem &mf2, const VEC2 &U2,
236  mesh_region rg = mesh_region::all_convexes()) {
237  return sqrt(asm_H1_semi_dist_sqr(mim, mf1, U1, mf2, U2, rg,
238  typename gmm::linalg_traits<VEC1>::value_type()));
239  }
240 
241  template<typename VEC1, typename VEC2, typename T>
242  inline scalar_type asm_H1_semi_dist_sqr
243  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
244  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, T) {
245  ga_workspace workspace;
246  model_real_plain_vector UU1(mf1.nb_dof()), UU2(mf2.nb_dof());
247  gmm::copy(U1, UU1); gmm::copy(U2, UU2);
248  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
249  workspace.add_fem_variable("u1", mf1, Iu1, UU1);
250  workspace.add_fem_variable("u2", mf2, Iu2, UU2);
251  workspace.add_expression("(Grad_u2-Grad_u1):(Grad_u2-Grad_u1)", mim, rg);
252  workspace.assembly(0);
253  return workspace.assembled_potential();
254  }
255 
256  inline scalar_type asm_H1_semi_dist_sqr
257  (const mesh_im &mim, const mesh_fem &mf1, const model_real_plain_vector &U1,
258  const mesh_fem &mf2, const model_real_plain_vector &U2,
259  mesh_region rg, scalar_type) {
260  ga_workspace workspace;
261  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
262  workspace.add_fem_variable("u1", mf1, Iu1, U1);
263  workspace.add_fem_variable("u2", mf2, Iu2, U2);
264  workspace.add_expression("(Grad_u2-Grad_u1):(Grad_u2-Grad_u1)", mim, rg);
265  workspace.assembly(0);
266  return workspace.assembled_potential();
267  }
268 
269  template<typename VEC1, typename VEC2, typename T>
270  inline scalar_type asm_H1_semi_dist_sqr
271  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
272  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, std::complex<T>) {
273  ga_workspace workspace;
274  model_real_plain_vector UU1R(mf1.nb_dof()), UU2R(mf2.nb_dof());
275  model_real_plain_vector UU1I(mf1.nb_dof()), UU2I(mf2.nb_dof());
276  gmm::copy(gmm::real_part(U1), UU1R); gmm::copy(gmm::imag_part(U1), UU1I);
277  gmm::copy(gmm::real_part(U2), UU2R); gmm::copy(gmm::imag_part(U2), UU2I);
278  gmm::sub_interval Iu1r(0, mf1.nb_dof()), Iu2r(mf1.nb_dof(), mf2.nb_dof());
279  gmm::sub_interval Iu1i(Iu2r.last(), mf1.nb_dof());
280  gmm::sub_interval Iu2i(Iu1i.last(), mf2.nb_dof());
281  workspace.add_fem_variable("u1", mf1, Iu1r, UU1R);
282  workspace.add_fem_variable("u2", mf2, Iu2r, UU2R);
283  workspace.add_fem_variable("v1", mf1, Iu1i, UU1I);
284  workspace.add_fem_variable("v2", mf2, Iu2i, UU2I);
285  workspace.add_expression("(Grad_u2-Grad_u1):(Grad_u2-Grad_u1)"
286  "+ (Grad_v2-Grad_v1):(Grad_v2-Grad_v1)", mim, rg);
287  workspace.assembly(0);
288  return workspace.assembled_potential();
289  }
290 
291  /**
292  compute the H1 norm of U.
293  @ingroup asm
294  */
295 
296  /**
297  compute @f$\|\nabla U\|_2@f$, U might be real or complex
298  @ingroup asm
299  */
300  template<typename VEC>
301  scalar_type asm_H1_norm
302  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
303  const mesh_region &rg = mesh_region::all_convexes()) {
304  typedef typename gmm::linalg_traits<VEC>::value_type T;
305  return sqrt(asm_H1_norm_sqr(mim, mf, U, rg, T()));
306  }
307 
308  template<typename VEC, typename T>
309  inline scalar_type asm_H1_norm_sqr
310  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
311  const mesh_region &rg_, T) {
312  ga_workspace workspace;
313  model_real_plain_vector UU(mf.nb_dof()); gmm::copy(U, UU);
314  gmm::sub_interval Iu(0, mf.nb_dof());
315  workspace.add_fem_variable("u", mf, Iu, UU);
316  workspace.add_expression("u.u + Grad_u:Grad_u", mim, rg_);
317  workspace.assembly(0);
318  return workspace.assembled_potential();
319  }
320 
321  inline scalar_type asm_H1_norm_sqr
322  (const mesh_im &mim, const mesh_fem &mf, const model_real_plain_vector &U,
323  const mesh_region &rg_, scalar_type) {
324  ga_workspace workspace;
325  gmm::sub_interval Iu(0, mf.nb_dof());
326  workspace.add_fem_variable("u", mf, Iu, U);
327  workspace.add_expression("u.u + Grad_u:Grad_u", mim, rg_);
328  workspace.assembly(0);
329  return workspace.assembled_potential();
330  }
331 
332 
333  template<typename VEC, typename T>
334  inline scalar_type asm_H1_norm_sqr
335  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
336  const mesh_region &rg, std::complex<T>) {
337  ga_workspace workspace;
338  model_real_plain_vector UUR(mf.nb_dof()), UUI(mf.nb_dof());
339  gmm::copy(gmm::real_part(U), UUR);
340  gmm::copy(gmm::imag_part(U), UUI);
341  gmm::sub_interval Iur(0, mf.nb_dof()), Iui(mf.nb_dof(), mf.nb_dof());
342  workspace.add_fem_variable("u", mf, Iur, UUR);
343  workspace.add_fem_variable("v", mf, Iui, UUI);
344  workspace.add_expression("u.u+v.v + Grad_u:Grad_u+Grad_v:Grad_v", mim, rg);
345  workspace.assembly(0);
346  return workspace.assembled_potential();
347  }
348 
349  /**
350  Compute the H1 distance between U1 and U2
351  @ingroup asm
352  */
353  template<typename VEC1, typename VEC2>
354  inline scalar_type asm_H1_dist
355  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
356  const mesh_fem &mf2, const VEC2 &U2,
357  mesh_region rg = mesh_region::all_convexes()) {
358  return sqrt(asm_H1_dist_sqr(mim, mf1, U1, mf2, U2, rg,
359  typename gmm::linalg_traits<VEC1>::value_type()));
360  }
361 
362  template<typename VEC1, typename VEC2, typename T>
363  inline scalar_type asm_H1_dist_sqr
364  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
365  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, T) {
366  ga_workspace workspace;
367  model_real_plain_vector UU1(mf1.nb_dof()), UU2(mf2.nb_dof());
368  gmm::copy(U1, UU1); gmm::copy(U2, UU2);
369  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
370  workspace.add_fem_variable("u1", mf1, Iu1, UU1);
371  workspace.add_fem_variable("u2", mf2, Iu2, UU2);
372  workspace.add_expression("(u2-u1).(u2-u1)"
373  "+ (Grad_u2-Grad_u1):(Grad_u2-Grad_u1)", mim, rg);
374  workspace.assembly(0);
375  return workspace.assembled_potential();
376  }
377 
378  inline scalar_type asm_H1_dist_sqr
379  (const mesh_im &mim, const mesh_fem &mf1, const model_real_plain_vector &U1,
380  const mesh_fem &mf2, const model_real_plain_vector &U2,
381  mesh_region rg, scalar_type) {
382  ga_workspace workspace;
383  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
384  workspace.add_fem_variable("u1", mf1, Iu1, U1);
385  workspace.add_fem_variable("u2", mf2, Iu2, U2);
386  workspace.add_expression("(u2-u1).(u2-u1)"
387  "+ (Grad_u2-Grad_u1):(Grad_u2-Grad_u1)", mim, rg);
388  workspace.assembly(0);
389  return workspace.assembled_potential();
390  }
391 
392  template<typename VEC1, typename VEC2, typename T>
393  inline scalar_type asm_H1_dist_sqr
394  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
395  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, std::complex<T>) {
396  ga_workspace workspace;
397  model_real_plain_vector UU1R(mf1.nb_dof()), UU2R(mf2.nb_dof());
398  model_real_plain_vector UU1I(mf1.nb_dof()), UU2I(mf2.nb_dof());
399  gmm::copy(gmm::real_part(U1), UU1R); gmm::copy(gmm::imag_part(U1), UU1I);
400  gmm::copy(gmm::real_part(U2), UU2R); gmm::copy(gmm::imag_part(U2), UU2I);
401  gmm::sub_interval Iu1r(0, mf1.nb_dof()), Iu2r(mf1.nb_dof(), mf2.nb_dof());
402  gmm::sub_interval Iu1i(Iu2r.last(), mf1.nb_dof());
403  gmm::sub_interval Iu2i(Iu1i.last(), mf2.nb_dof());
404  workspace.add_fem_variable("u1", mf1, Iu1r, UU1R);
405  workspace.add_fem_variable("u2", mf2, Iu2r, UU2R);
406  workspace.add_fem_variable("v1", mf1, Iu1i, UU1I);
407  workspace.add_fem_variable("v2", mf2, Iu2i, UU2I);
408  workspace.add_expression("(u2-u1).(u2-u1) + (v2-v1).(v2-v1)"
409  "+ (Grad_u2-Grad_u1):(Grad_u2-Grad_u1)"
410  "+ (Grad_v2-Grad_v1):(Grad_v2-Grad_v1)", mim, rg);
411  workspace.assembly(0);
412  return workspace.assembled_potential();
413  }
414 
415  /**
416  compute @f$\|Hess U\|_2@f$, U might be real or complex. For C^1 elements
417  @ingroup asm
418  */
419 
420  template<typename VEC>
421  scalar_type asm_H2_semi_norm
422  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
423  const mesh_region &rg = mesh_region::all_convexes()) {
424  typedef typename gmm::linalg_traits<VEC>::value_type T;
425  return sqrt(asm_H2_semi_norm_sqr(mim, mf, U, rg, T()));
426  }
427 
428  template<typename VEC, typename T>
429  inline scalar_type asm_H2_semi_norm_sqr
430  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
431  const mesh_region &rg_, T) {
432  ga_workspace workspace;
433  model_real_plain_vector UU(mf.nb_dof()); gmm::copy(U, UU);
434  gmm::sub_interval Iu(0, mf.nb_dof());
435  workspace.add_fem_variable("u", mf, Iu, UU);
436  workspace.add_expression("Hess_u:Hess_u", mim, rg_);
437  workspace.assembly(0);
438  return workspace.assembled_potential();
439  }
440 
441  inline scalar_type asm_H2_semi_norm_sqr
442  (const mesh_im &mim, const mesh_fem &mf, const model_real_plain_vector &U,
443  const mesh_region &rg_, scalar_type) {
444  ga_workspace workspace;
445  gmm::sub_interval Iu(0, mf.nb_dof());
446  workspace.add_fem_variable("u", mf, Iu, U);
447  workspace.add_expression("Hess_u:Hess_u", mim, rg_);
448  workspace.assembly(0);
449  return workspace.assembled_potential();
450  }
451 
452 
453  template<typename VEC, typename T>
454  inline scalar_type asm_H2_semi_norm_sqr
455  (const mesh_im &mim, const mesh_fem &mf, const VEC &U,
456  const mesh_region &rg, std::complex<T>) {
457  ga_workspace workspace;
458  model_real_plain_vector UUR(mf.nb_dof()), UUI(mf.nb_dof());
459  gmm::copy(gmm::real_part(U), UUR);
460  gmm::copy(gmm::imag_part(U), UUI);
461  gmm::sub_interval Iur(0, mf.nb_dof()), Iui(mf.nb_dof(), mf.nb_dof());
462  workspace.add_fem_variable("u", mf, Iur, UUR);
463  workspace.add_fem_variable("v", mf, Iui, UUI);
464  workspace.add_expression("Hess_u:Hess_u + Hess_v:Hess_v", mim, rg);
465  workspace.assembly(0);
466  return workspace.assembled_potential();
467  }
468 
469 
470  template<typename VEC1, typename VEC2>
471  inline scalar_type asm_H2_semi_dist
472  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
473  const mesh_fem &mf2, const VEC2 &U2,
474  mesh_region rg = mesh_region::all_convexes()) {
475  return sqrt(asm_H2_semi_dist_sqr(mim, mf1, U1, mf2, U2, rg,
476  typename gmm::linalg_traits<VEC1>::value_type()));
477  }
478 
479  template<typename VEC1, typename VEC2, typename T>
480  inline scalar_type asm_H2_semi_dist_sqr
481  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
482  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, T) {
483  ga_workspace workspace;
484  model_real_plain_vector UU1(mf1.nb_dof()), UU2(mf2.nb_dof());
485  gmm::copy(U1, UU1); gmm::copy(U2, UU2);
486  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
487  workspace.add_fem_variable("u1", mf1, Iu1, UU1);
488  workspace.add_fem_variable("u2", mf2, Iu2, UU2);
489  workspace.add_expression("(Hess_u2-Hess_u1):(Hess_u2-Hess_u1)", mim, rg);
490  workspace.assembly(0);
491  return workspace.assembled_potential();
492  }
493 
494  inline scalar_type asm_H2_semi_dist_sqr
495  (const mesh_im &mim, const mesh_fem &mf1, const model_real_plain_vector &U1,
496  const mesh_fem &mf2, const model_real_plain_vector &U2,
497  mesh_region rg, scalar_type) {
498  ga_workspace workspace;
499  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(mf1.nb_dof(), mf2.nb_dof());
500  workspace.add_fem_variable("u1", mf1, Iu1, U1);
501  workspace.add_fem_variable("u2", mf2, Iu2, U2);
502  workspace.add_expression("(Hess_u2-Hess_u1):(Hess_u2-Hess_u1)", mim, rg);
503  workspace.assembly(0);
504  return workspace.assembled_potential();
505  }
506 
507  template<typename VEC1, typename VEC2, typename T>
508  inline scalar_type asm_H2_semi_dist_sqr
509  (const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1,
510  const mesh_fem &mf2, const VEC2 &U2, mesh_region rg, std::complex<T>) {
511  ga_workspace workspace;
512  model_real_plain_vector UU1R(mf1.nb_dof()), UU2R(mf2.nb_dof());
513  model_real_plain_vector UU1I(mf1.nb_dof()), UU2I(mf2.nb_dof());
514  gmm::copy(gmm::real_part(U1), UU1R); gmm::copy(gmm::imag_part(U1), UU1I);
515  gmm::copy(gmm::real_part(U2), UU2R); gmm::copy(gmm::imag_part(U2), UU2I);
516  gmm::sub_interval Iu1r(0, mf1.nb_dof()), Iu2r(mf1.nb_dof(), mf2.nb_dof());
517  gmm::sub_interval Iu1i(Iu2r.last(), mf1.nb_dof());
518  gmm::sub_interval Iu2i(Iu1i.last(), mf2.nb_dof());
519  workspace.add_fem_variable("u1", mf1, Iu1r, UU1R);
520  workspace.add_fem_variable("u2", mf2, Iu2r, UU2R);
521  workspace.add_fem_variable("v1", mf1, Iu1i, UU1I);
522  workspace.add_fem_variable("v2", mf2, Iu2i, UU2I);
523  workspace.add_expression("(Hess_u2-Hess_u1):(Hess_u2-Hess_u1)"
524  "+ (Hess_v2-Hess_v1):(Hess_v2-Hess_v1)", mim, rg);
525  workspace.assembly(0);
526  return workspace.assembled_potential();
527  }
528 
529  /**
530  compute the H2 norm of U (for C^1 elements).
531  @ingroup asm
532  */
533  template<typename VEC>
534  scalar_type asm_H2_norm(const mesh_im &mim, const mesh_fem &mf,
535  const VEC &U,
536  const mesh_region &rg
537  = mesh_region::all_convexes()) {
538  typedef typename gmm::linalg_traits<VEC>::value_type T;
539  return sqrt(asm_H1_norm_sqr(mim, mf, U, rg, T())
540  + asm_H2_semi_norm_sqr(mim, mf, U, rg, T()));
541  }
542 
543  /**
544  Compute the H2 distance between U1 and U2
545  @ingroup asm
546  */
547  template<typename VEC1, typename VEC2>
548  scalar_type asm_H2_dist(const mesh_im &mim,
549  const mesh_fem &mf1, const VEC1 &U1,
550  const mesh_fem &mf2, const VEC2 &U2,
551  const mesh_region &rg
552  = mesh_region::all_convexes()) {
553  typedef typename gmm::linalg_traits<VEC1>::value_type T;
554  return sqrt(asm_H1_dist_sqr(mim,mf1,U1,mf2,U2,rg,T()) +
555  asm_H2_semi_dist_sqr(mim,mf1,U1,mf2,U2,rg,T()));
556  }
557 
558  /*
559  assembly of a matrix with 1 parameter (real or complex)
560  (the most common here for the assembly routines below)
561  */
562  template <typename MAT, typename VECT>
563  inline void asm_real_or_complex_1_param_mat
564  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem *mf_data,
565  const VECT &A, const mesh_region &rg, const char *assembly_description) {
566  asm_real_or_complex_1_param_mat_
567  (M, mim, mf_u, mf_data, A, rg, assembly_description,
568  typename gmm::linalg_traits<VECT>::value_type());
569  }
570 
571  /* real version */
572  template<typename MAT, typename VECT, typename T>
573  inline void asm_real_or_complex_1_param_mat_
574  (const MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
575  const mesh_fem *mf_data, const VECT &A, const mesh_region &rg,
576  const char *assembly_description, T) {
577  ga_workspace workspace;
578  gmm::sub_interval Iu(0, mf_u.nb_dof());
579  base_vector u(mf_u.nb_dof()), AA(gmm::vect_size(A));
580  gmm::copy(A, AA);
581  workspace.add_fem_variable("u", mf_u, Iu, u);
582  if (mf_data)
583  workspace.add_fem_constant("A", *mf_data, AA);
584  else
585  workspace.add_fixed_size_constant("A", AA);
586  workspace.add_expression(assembly_description, mim, rg);
587  workspace.assembly(2);
588  if (gmm::mat_nrows(workspace.assembled_matrix()))
589  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
590  }
591 
592  inline void asm_real_or_complex_1_param_mat_
593  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf_u,
594  const mesh_fem *mf_data, const model_real_plain_vector &A,
595  const mesh_region &rg,
596  const char *assembly_description, scalar_type) {
597  ga_workspace workspace;
598  gmm::sub_interval Iu(0, mf_u.nb_dof());
599  base_vector u(mf_u.nb_dof());
600  workspace.add_fem_variable("u", mf_u, Iu, u);
601  if (mf_data)
602  workspace.add_fem_constant("A", *mf_data, A);
603  else
604  workspace.add_fixed_size_constant("A", A);
605  workspace.add_expression(assembly_description, mim, rg);
606  workspace.set_assembled_matrix(M);
607  workspace.assembly(2);
608  }
609 
610  /* complex version */
611  template<typename MAT, typename VECT, typename T>
612  inline void asm_real_or_complex_1_param_mat_
613  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem *mf_data,
614  const VECT &A, const mesh_region &rg, const char *assembly_description,
615  std::complex<T>) {
616  asm_real_or_complex_1_param_mat_(gmm::real_part(M),mim,mf_u,mf_data,
617  gmm::real_part(A),rg,
618  assembly_description, T());
619  asm_real_or_complex_1_param_mat_(gmm::imag_part(M),mim,mf_u,mf_data,
620  gmm::imag_part(A),rg,
621  assembly_description, T());
622  }
623 
624  /*
625  assembly of a vector with 1 parameter (real or complex)
626  (the most common here for the assembly routines below)
627  */
628  template <typename MAT, typename VECT>
629  inline void asm_real_or_complex_1_param_vec
630  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem *mf_data,
631  const VECT &A, const mesh_region &rg, const char *assembly_description) {
632  asm_real_or_complex_1_param_vec_
633  (M, mim, mf_u, mf_data, A, rg, assembly_description,
634  typename gmm::linalg_traits<VECT>::value_type());
635  }
636 
637  /* real version */
638  template<typename VECTA, typename VECT, typename T>
639  inline void asm_real_or_complex_1_param_vec_
640  (const VECTA &V, const mesh_im &mim, const mesh_fem &mf_u,
641  const mesh_fem *mf_data, const VECT &A, const mesh_region &rg,
642  const char *assembly_description, T) {
643  ga_workspace workspace;
644  gmm::sub_interval Iu(0, mf_u.nb_dof());
645  base_vector u(mf_u.nb_dof()), AA(gmm::vect_size(A));
646  gmm::copy(A, AA);
647  workspace.add_fem_variable("u", mf_u, Iu, u);
648  if (mf_data)
649  workspace.add_fem_constant("A", *mf_data, AA);
650  else
651  workspace.add_fixed_size_constant("A", AA);
652  workspace.add_expression(assembly_description, mim, rg);
653  workspace.assembly(1);
654  if (gmm::vect_size(workspace.assembled_vector()))
655  gmm::add(workspace.assembled_vector(), const_cast<VECTA &>(V));
656  }
657 
658  inline void asm_real_or_complex_1_param_vec_
659  (model_real_plain_vector &V, const mesh_im &mim, const mesh_fem &mf_u,
660  const mesh_fem *mf_data, const model_real_plain_vector &A,
661  const mesh_region &rg,
662  const char *assembly_description, scalar_type) {
663  ga_workspace workspace;
664  gmm::sub_interval Iu(0, mf_u.nb_dof());
665  base_vector u(mf_u.nb_dof());
666  workspace.add_fem_variable("u", mf_u, Iu, u);
667  if (mf_data)
668  workspace.add_fem_constant("A", *mf_data, A);
669  else
670  workspace.add_fixed_size_constant("A", A);
671  workspace.add_expression(assembly_description, mim, rg);
672  workspace.set_assembled_vector(V);
673  workspace.assembly(1);
674  }
675 
676  /* complex version */
677  template<typename MAT, typename VECT, typename T>
678  inline void asm_real_or_complex_1_param_vec_
679  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem *mf_data,
680  const VECT &A, const mesh_region &rg,const char *assembly_description,
681  std::complex<T>) {
682  asm_real_or_complex_1_param_vec_(gmm::real_part(M),mim,mf_u,mf_data,
683  gmm::real_part(A),rg,
684  assembly_description, T());
685  asm_real_or_complex_1_param_vec_(gmm::imag_part(M),mim,mf_u,mf_data,
686  gmm::imag_part(A),rg,
687  assembly_description, T());
688  }
689 
690  /**
691  generic mass matrix assembly (on the whole mesh or on the specified
692  convex set or boundary)
693  @ingroup asm
694  */
695  template<typename MAT>
696  inline void asm_mass_matrix
697  (const MAT &M, const mesh_im &mim, const mesh_fem &mf1,
698  const mesh_region &rg = mesh_region::all_convexes()) {
699 
700  ga_workspace workspace;
701  gmm::sub_interval Iu1(0, mf1.nb_dof());
702  base_vector u1(mf1.nb_dof());
703  workspace.add_fem_variable("u1", mf1, Iu1, u1);
704  workspace.add_expression("Test_u1:Test2_u1", mim, rg);
705  workspace.assembly(2);
706  if (gmm::mat_nrows(workspace.assembled_matrix()))
707  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
708  }
709 
710  inline void asm_mass_matrix
711  (model_real_sparse_matrix &M, const mesh_im &mim,
712  const mesh_fem &mf1,
713  const mesh_region &rg = mesh_region::all_convexes()) {
714  ga_workspace workspace;
715  gmm::sub_interval Iu1(0, mf1.nb_dof());
716  base_vector u1(mf1.nb_dof());
717  workspace.add_fem_variable("u1", mf1, Iu1, u1);
718  workspace.add_expression("Test_u1:Test2_u1", mim, rg);
719  workspace.set_assembled_matrix(M);
720  workspace.assembly(2);
721  }
722 
723  /**
724  * generic mass matrix assembly (on the whole mesh or on the specified
725  * boundary)
726  */
727 
728  template<typename MAT>
729  inline void asm_mass_matrix
730  (const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_fem &mf2,
731  const mesh_region &rg = mesh_region::all_convexes()) {
732  ga_workspace workspace;
733  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(Iu1.last(), mf2.nb_dof());
734  base_vector u1(mf1.nb_dof()), u2(mf2.nb_dof());
735  workspace.add_fem_variable("u1", mf1, Iu1, u1);
736  workspace.add_fem_variable("u2", mf2, Iu2, u2);
737  workspace.add_expression("Test_u1:Test2_u2", mim, rg);
738  workspace.assembly(2);
739  if (gmm::mat_nrows(workspace.assembled_matrix()))
740  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(), Iu1, Iu2),
741  const_cast<MAT &>(M));
742  }
743 
744  inline void asm_mass_matrix
745  (model_real_sparse_matrix &M, const mesh_im &mim,
746  const mesh_fem &mf1, const mesh_fem &mf2,
747  const mesh_region &rg = mesh_region::all_convexes()) {
748  ga_workspace workspace;
749  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(0, mf2.nb_dof());
750  base_vector u1(mf1.nb_dof()), u2(mf2.nb_dof());
751  workspace.add_fem_variable("u1", mf1, Iu1, u1);
752  workspace.add_fem_variable("u2", mf2, Iu2, u2);
753  workspace.add_expression("Test_u1:Test2_u2", mim, rg);
754  workspace.set_assembled_matrix(M);
755  workspace.assembly(2);
756  }
757 
758  /**
759  generic mass matrix assembly with an additional parameter
760  (on the whole mesh or on the specified boundary)
761  @ingroup asm
762  */
763  template<typename MAT, typename VECT>
764  inline void asm_mass_matrix_param
765  (const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_fem &mf2,
766  const mesh_fem &mf_data, const VECT &A,
767  const mesh_region &rg = mesh_region::all_convexes()) {
768 
769  ga_workspace workspace;
770  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(Iu1.last(), mf2.nb_dof());
771  base_vector u1(mf1.nb_dof()), u2(mf2.nb_dof()), AA(mf_data.nb_dof());
772  gmm::copy(A, AA);
773  workspace.add_fem_variable("u1", mf1, Iu1, u1);
774  workspace.add_fem_variable("u2", mf2, Iu2, u2);
775  workspace.add_fem_constant("A", mf_data, AA);
776  workspace.add_expression("(A*Test_u1).Test2_u2", mim, rg);
777  workspace.assembly(2);
778  if (gmm::mat_nrows(workspace.assembled_matrix()))
779  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(), Iu1, Iu2),
780  const_cast<MAT &>(M));
781  }
782 
783  inline void asm_mass_matrix_param
784  (model_real_sparse_matrix &M, const mesh_im &mim,
785  const mesh_fem &mf1, const mesh_fem &mf2,
786  const mesh_fem &mf_data, const model_real_plain_vector &A,
787  const mesh_region &rg = mesh_region::all_convexes()) {
788 
789  ga_workspace workspace;
790  gmm::sub_interval Iu1(0, mf1.nb_dof()), Iu2(0, mf2.nb_dof());
791  base_vector u1(mf1.nb_dof()), u2(mf2.nb_dof());
792  workspace.add_fem_variable("u1", mf1, Iu1, u1);
793  workspace.add_fem_variable("u2", mf2, Iu2, u2);
794  workspace.add_fem_constant("A", mf_data, A);
795  workspace.add_expression("(A*Test_u1).Test2_u2", mim, rg);
796  workspace.set_assembled_matrix(M);
797  workspace.assembly(2);
798  }
799 
800 
801 
802  /**
803  generic mass matrix assembly with an additional parameter
804  (on the whole mesh or on the specified boundary)
805  @ingroup asm
806  */
807  template<typename MAT, typename VECT>
808  void asm_mass_matrix_param
809  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_data,
810  const VECT &F, const mesh_region &rg = mesh_region::all_convexes()) {
811  asm_real_or_complex_1_param_mat
812  (M, mim, mf_u, &mf_data, F, rg, "(A*Test_u):Test2_u");
813  }
814 
815  /**
816  generic mass matrix assembly with an additional constant parameter
817  (on the whole mesh or on the specified boundary)
818  @ingroup asm
819  */
820  template<typename MAT, typename VECT>
821  void asm_mass_matrix_homogeneous_param
822  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
823  const VECT &F, const mesh_region &rg = mesh_region::all_convexes()) {
824  asm_real_or_complex_1_param_mat
825  (M, mim, mf_u, 0, F, rg, "(A*Test_u):Test2_u");
826  }
827 
828  /**
829  source term (for both volumic sources and boundary (Neumann) sources).
830  @ingroup asm
831  */
832  template<typename VECT1, typename VECT2>
833  void asm_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf,
834  const mesh_fem &mf_data, const VECT2 &F,
835  const mesh_region &rg = mesh_region::all_convexes()) {
836  GMM_ASSERT1(mf_data.get_qdim() == 1 ||
837  mf_data.get_qdim() == mf.get_qdim(),
838  "invalid data mesh fem (same Qdim or Qdim=1 required)");
839  asm_real_or_complex_1_param_vec
840  (const_cast<VECT1 &>(B), mim, mf, &mf_data, F, rg, "A:Test_u");
841  }
842 
843  /**
844  source term (for both volumic sources and boundary (Neumann) sources).
845  For an homogeneous term.
846  @ingroup asm
847  */
848  template<typename VECT1, typename VECT2>
849  void asm_homogeneous_source_term(const VECT1 &B, const mesh_im &mim,
850  const mesh_fem &mf, const VECT2 &F,
851  const mesh_region &rg = mesh_region::all_convexes()) {
852  asm_real_or_complex_1_param_vec
853  (const_cast<VECT1 &>(B), mim, mf, 0, F, rg, "A:Test_u");
854  }
855 
856  /**
857  Normal source term (for boundary (Neumann) condition).
858  @ingroup asm
859  */
860  template<typename VECT1, typename VECT2>
861  void asm_normal_source_term(VECT1 &B, const mesh_im &mim,
862  const mesh_fem &mf,
863  const mesh_fem &mf_data, const VECT2 &F,
864  const mesh_region &rg) {
865  asm_real_or_complex_1_param_vec(B, mim, mf, &mf_data, F, rg,
866  "(Reshape(A, qdim(u), meshdim).Normal):Test_u");
867  }
868 
869  /**
870  Homogeneous normal source term (for boundary (Neumann) condition).
871  @ingroup asm
872  */
873  template<typename VECT1, typename VECT2>
874  void asm_homogeneous_normal_source_term
875  (VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const VECT2 &F,
876  const mesh_region &rg)
877  { asm_real_or_complex_1_param_vec(B, mim, mf, 0,F,rg,
878  "(Reshape(A, qdim(u), meshdim).Normal):Test_u"); }
879 
880  /**
881  assembly of @f$\int{qu.v}@f$
882 
883  (if @f$u@f$ is a vector field of size @f$N@f$, @f$q@f$ is a square
884  matrix @f$N\times N@f$ used by assem_general_boundary_conditions
885 
886  convention: Q is of the form
887  Q1_11 Q2_11 ..... Qn_11
888  Q1_21 Q2_21 ..... Qn_21
889  Q1_12 Q2_12 ..... Qn_12
890  Q1_22 Q2_22 ..... Qn_22
891  if N = 2, and mf_d has n/N degree of freedom
892 
893  Q is a vector, so the matrix is assumed to be stored by columns
894  (fortran style)
895 
896  Works for both volumic assembly and boundary assembly
897  @ingroup asm
898  */
899  template<typename MAT, typename VECT>
900  void asm_qu_term(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
901  const mesh_fem &mf_d, const VECT &Q,
902  const mesh_region &rg) {
903  if (mf_d.get_qdim() == 1 && gmm::vect_size(Q) > mf_d.nb_dof())
904  asm_real_or_complex_1_param_mat
905  (M, mim,mf_u,&mf_d,Q,rg,"(Reshape(A,qdim(u),qdim(u)).Test_u):Test2_u");
906  else if (mf_d.get_qdim() == mf_u.get_qdim())
907  asm_real_or_complex_1_param_mat
908  (M, mim, mf_u, &mf_d, Q, rg, "(A*Test_u):Test2_u");
909  else GMM_ASSERT1(false, "invalid data mesh fem");
910  }
911 
912  template<typename MAT, typename VECT>
913  void asm_homogeneous_qu_term(MAT &M, const mesh_im &mim,
914  const mesh_fem &mf_u, const VECT &Q,
915  const mesh_region &rg) {
916  if (gmm::vect_size(Q) == 1)
917  asm_real_or_complex_1_param_mat
918  (M, mim,mf_u,0,Q,rg,"(A*Test_u):Test2_u");
919  else
920  asm_real_or_complex_1_param_mat
921  (M, mim,mf_u,0,Q,rg,"(Reshape(A,qdim(u),qdim(u)).Test_u):Test2_u");
922  }
923 
924  /**
925  Stiffness matrix for linear elasticity, with Lamé coefficients
926  @ingroup asm
927  */
928  template<class MAT, class VECT>
929  inline void asm_stiffness_matrix_for_linear_elasticity
930  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
931  const mesh_fem &mf_data, const VECT &LAMBDA, const VECT &MU,
932  const mesh_region &rg = mesh_region::all_convexes()) {
933  ga_workspace workspace;
934  gmm::sub_interval Iu(0, mf.nb_dof());
935  base_vector u(mf.nb_dof()), lambda(gmm::vect_size(LAMBDA));
936  base_vector mu(gmm::vect_size(MU));
937  gmm::copy(LAMBDA, lambda); gmm::copy(MU, mu);
938  workspace.add_fem_variable("u", mf, Iu, u);
939  workspace.add_fem_constant("lambda", mf_data, lambda);
940  workspace.add_fem_constant("mu", mf_data, mu);
941  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
942  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
943  workspace.assembly(2);
944  if (gmm::mat_nrows(workspace.assembled_matrix()))
945  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
946  }
947 
948  inline void asm_stiffness_matrix_for_linear_elasticity
949  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf,
950  const mesh_fem &mf_data, const model_real_plain_vector &LAMBDA,
951  const model_real_plain_vector &MU,
952  const mesh_region &rg = mesh_region::all_convexes()) {
953  ga_workspace workspace;
954  gmm::sub_interval Iu(0, mf.nb_dof());
955  base_vector u(mf.nb_dof());
956  workspace.add_fem_variable("u", mf, Iu, u);
957  workspace.add_fem_constant("lambda", mf_data, LAMBDA);
958  workspace.add_fem_constant("mu", mf_data, MU);
959  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
960  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
961  workspace.set_assembled_matrix(M);
962  workspace.assembly(2);
963  }
964 
965 
966  /**
967  Stiffness matrix for linear elasticity, with constant Lamé coefficients
968  @ingroup asm
969  */
970  template<class MAT, class VECT>
971  inline void asm_stiffness_matrix_for_homogeneous_linear_elasticity
972  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
973  const VECT &LAMBDA, const VECT &MU,
974  const mesh_region &rg = mesh_region::all_convexes()) {
975  ga_workspace workspace;
976  gmm::sub_interval Iu(0, mf.nb_dof());
977  base_vector u(mf.nb_dof()), lambda(gmm::vect_size(LAMBDA));
978  base_vector mu(gmm::vect_size(MU));
979  gmm::copy(LAMBDA, lambda); gmm::copy(MU, mu);
980  workspace.add_fem_variable("u", mf, Iu, u);
981  workspace.add_fixed_size_constant("lambda", lambda);
982  workspace.add_fixed_size_constant("mu", mu);
983  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
984  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
985  workspace.assembly(2);
986  if (gmm::mat_nrows(workspace.assembled_matrix()))
987  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
988  }
989 
990 
991  inline void asm_stiffness_matrix_for_homogeneous_linear_elasticity
992  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf,
993  const model_real_plain_vector &LAMBDA, const model_real_plain_vector &MU,
994  const mesh_region &rg = mesh_region::all_convexes()) {
995  ga_workspace workspace;
996  gmm::sub_interval Iu(0, mf.nb_dof());
997  base_vector u(mf.nb_dof());
998  workspace.add_fem_variable("u", mf, Iu, u);
999  workspace.add_fixed_size_constant("lambda", LAMBDA);
1000  workspace.add_fixed_size_constant("mu", MU);
1001  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
1002  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
1003  workspace.set_assembled_matrix(M);
1004  workspace.assembly(2);
1005  }
1006 
1007  /**
1008  Stiffness matrix for linear elasticity, with a general Hooke
1009  tensor. This is more a demonstration of generic assembly than
1010  something useful !
1011 
1012  Note that this function is just an alias for
1013  asm_stiffness_matrix_for_vector_elliptic.
1014 
1015  @ingroup asm
1016  */
1017  template<typename MAT, typename VECT> void
1018  asm_stiffness_matrix_for_linear_elasticity_Hooke
1019  (MAT &RM, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1020  const VECT &H, const mesh_region &rg = mesh_region::all_convexes()) {
1021  asm_stiffness_matrix_for_vector_elliptic(RM, mim, mf, mf_data, H, rg);
1022  }
1023 
1024  /**
1025  Build the mixed pressure term @f$ B = - \int p.div u @f$
1026 
1027  @ingroup asm
1028  */
1029  template<typename MAT>
1030  inline void asm_stokes_B(const MAT &B, const mesh_im &mim,
1031  const mesh_fem &mf_u, const mesh_fem &mf_p,
1032  const mesh_region &rg=mesh_region::all_convexes()) {
1033  ga_workspace workspace;
1034  gmm::sub_interval Iu(0, mf_u.nb_dof()), Ip(Iu.last(), mf_p.nb_dof());
1035  base_vector u(mf_u.nb_dof()), p(mf_p.nb_dof());
1036  workspace.add_fem_variable("u", mf_u, Iu, u);
1037  workspace.add_fem_variable("p", mf_p, Ip, p);
1038  workspace.add_expression("Test_p*Div_Test2_u", mim, rg);
1039  workspace.assembly(2);
1040  if (gmm::mat_nrows(workspace.assembled_matrix()))
1041  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(), Ip, Iu),
1042  const_cast<MAT &>(B));
1043  }
1044 
1045  inline void asm_stokes_B(model_real_sparse_matrix &B, const mesh_im &mim,
1046  const mesh_fem &mf_u, const mesh_fem &mf_p,
1047  const mesh_region &rg=mesh_region::all_convexes()) {
1048  ga_workspace workspace;
1049  gmm::sub_interval Iu(0, mf_u.nb_dof()), Ip(0, mf_p.nb_dof());
1050  base_vector u(mf_u.nb_dof()), p(mf_p.nb_dof());
1051  workspace.add_fem_variable("u", mf_u, Iu, u);
1052  workspace.add_fem_variable("p", mf_p, Ip, p);
1053  workspace.add_expression("Test_p*Div_Test2_u", mim, rg);
1054  workspace.set_assembled_matrix(B);
1055  workspace.assembly(2);
1056  }
1057 
1058 
1059  /**
1060  assembly of @f$\int_\Omega \nabla u.\nabla v@f$.
1061 
1062  @ingroup asm
1063  */
1064  template<typename MAT>
1065  inline void asm_stiffness_matrix_for_homogeneous_laplacian
1066  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
1067  const mesh_region &rg = mesh_region::all_convexes()) {
1068  ga_workspace workspace;
1069  gmm::sub_interval Iu(0, mf.nb_dof());
1070  base_vector u(mf.nb_dof());
1071  workspace.add_fem_variable("u", mf, Iu, u);
1072  workspace.add_expression("Grad_Test_u:Grad_Test2_u", mim, rg);
1073  workspace.assembly(2);
1074  if (gmm::mat_nrows(workspace.assembled_matrix()))
1075  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
1076  }
1077 
1078  inline void asm_stiffness_matrix_for_homogeneous_laplacian
1079  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf,
1080  const mesh_region &rg = mesh_region::all_convexes()) {
1081  ga_workspace workspace;
1082  gmm::sub_interval Iu(0, mf.nb_dof());
1083  base_vector u(mf.nb_dof());
1084  workspace.add_fem_variable("u", mf, Iu, u);
1085  workspace.add_expression("Grad_Test_u:Grad_Test2_u", mim, rg);
1086  workspace.set_assembled_matrix(M);
1087  workspace.assembly(2);
1088  }
1089 
1090  /**
1091  assembly of @f$\int_\Omega \nabla u.\nabla v@f$.
1092  @ingroup asm
1093  */
1094  template<typename MAT>
1095  inline void asm_stiffness_matrix_for_homogeneous_laplacian_componentwise
1096  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
1097  const mesh_region &rg = mesh_region::all_convexes()) {
1098  return asm_stiffness_matrix_for_homogeneous_laplacian(M, mim, mf, rg);
1099  }
1100 
1101  /**
1102  assembly of @f$\int_\Omega a(x)\nabla u.\nabla v@f$ , where @f$a(x)@f$
1103  is scalar.
1104  @ingroup asm
1105  */
1106  template<typename MAT, typename VECT>
1107  inline void asm_stiffness_matrix_for_laplacian
1108  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1109  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1110  GMM_ASSERT1(mf_data.get_qdim() == 1
1111  && gmm::vect_size(A) == mf_data.nb_dof(), "invalid data");
1112  asm_real_or_complex_1_param_mat
1113  (M, mim, mf, &mf_data, A, rg, "(A*Grad_Test_u):Grad_Test2_u");
1114  }
1115 
1116  /** The same as getfem::asm_stiffness_matrix_for_laplacian , but on
1117  each component of mf when mf has a qdim > 1
1118  */
1119  template<typename MAT, typename VECT>
1120  inline void asm_stiffness_matrix_for_laplacian_componentwise
1121  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1122  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1123  asm_stiffness_matrix_for_laplacian(M, mim, mf, mf_data, A, rg);
1124  }
1125 
1126  /**
1127  assembly of @f$\int_\Omega A(x)\nabla u.\nabla v@f$, where @f$A(x)@f$
1128  is a (symmetric positive definite) NxN matrix.
1129  Arguments:
1130  @param M a sparse matrix of dimensions mf.nb_dof() x mf.nb_dof()
1131 
1132  @param mim the mesh_im.
1133 
1134  @param mf : the mesh_fem that describes the solution, with
1135  @c mf.get_qdim() == @c N.
1136 
1137  @param mf_data the mesh_fem that describes the coefficients of @c A
1138  (@c mf_data.get_qdim() == 1).
1139 
1140  @param A a (very large) vector, which is a flattened (n x n x
1141  mf_data.nb_dof()) 3D array. For each dof of mf_data, it contains
1142  the n x n coefficients of @f$A@f$. As usual, the order is the
1143  "fortran-order", i.e. @c A = [A_11(dof1) A_21(dof1) A_31(dof1)
1144  A_12(dof1) A_22(dof1) ... A_33(dof) A_11(dof2)
1145  .... A_33(lastdof)]
1146 
1147  @ingroup asm
1148  */
1149  template<typename MAT, typename VECT>
1150  void asm_stiffness_matrix_for_scalar_elliptic
1151  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1152  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1153  asm_real_or_complex_1_param_mat
1154  (M, mim, mf, &mf_data, A, rg,
1155  "(Reshape(A,meshdim,meshdim)*Grad_Test_u):Grad_Test2_u");
1156  }
1157 
1158  /** The same but with a constant matrix
1159  */
1160  template<typename MAT, typename VECT>
1161  void asm_stiffness_matrix_for_homogeneous_scalar_elliptic
1162  (MAT &M, const mesh_im &mim, const mesh_fem &mf,
1163  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1164  asm_real_or_complex_1_param_mat
1165  (M, mim, mf, 0, A, rg,
1166  "(Reshape(A,meshdim,meshdim)*Grad_Test_u):Grad_Test2_u");
1167  }
1168 
1169  /** The same but on each component of mf when mf has a qdim > 1
1170  */
1171  template<typename MAT, typename VECT>
1172  void asm_stiffness_matrix_for_scalar_elliptic_componentwise
1173  (MAT &M, const mesh_im &mim, const mesh_fem &mf,
1174  const mesh_fem &mf_data, const VECT &A,
1175  const mesh_region &rg = mesh_region::all_convexes()) {
1176  asm_real_or_complex_1_param_mat
1177  (M, mim, mf, &mf_data, A, rg,
1178  "(Grad_Test_u*(Reshape(A,meshdim,meshdim)')):Grad_Test2_u");
1179  }
1180 
1181  /** The same but with a constant matrix
1182  */
1183  template<typename MAT, typename VECT>
1184  void asm_stiffness_matrix_for_homogeneous_scalar_elliptic_componentwise
1185  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A,
1186  const mesh_region &rg = mesh_region::all_convexes()) {
1187  asm_real_or_complex_1_param_mat
1188  (M, mim, mf, 0, A, rg,
1189  "(Grad_Test_u*(Reshape(A,meshdim,meshdim)')):Grad_Test2_u");
1190  }
1191 
1192 
1193  /**
1194  Assembly of @f$\int_\Omega A(x)\nabla u.\nabla v@f$, where @f$A(x)@f$
1195  is a NxNxQxQ (symmetric positive definite) tensor defined on mf_data.
1196  */
1197  template<typename MAT, typename VECT> void
1198  asm_stiffness_matrix_for_vector_elliptic
1199  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1200  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1201  /* M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v) */
1202  asm_real_or_complex_1_param_mat
1203  (M, mim, mf, &mf_data, A, rg,
1204  "(Reshape(A,qdim(u),meshdim,qdim(u),meshdim):Grad_Test_u):Grad_Test2_u");
1205  }
1206 
1207  /**
1208  Assembly of @f$\int_\Omega A(x)\nabla u.\nabla v@f$, where @f$A(x)@f$
1209  is a NxNxQxQ (symmetric positive definite) constant tensor.
1210  */
1211  template<typename MAT, typename VECT> void
1212  asm_stiffness_matrix_for_homogeneous_vector_elliptic
1213  (MAT &M, const mesh_im &mim, const mesh_fem &mf,
1214  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1215  /* M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v) */
1216  asm_real_or_complex_1_param_mat
1217  (M, mim, mf, 0, A, rg,
1218  "(Reshape(A,qdim(u),meshdim,qdim(u),meshdim):Grad_Test_u):Grad_Test2_u");
1219  }
1220 
1221  /**
1222  assembly of the term @f$\int_\Omega Kuv - \nabla u.\nabla v@f$,
1223  for the helmholtz equation (@f$\Delta u + k^2u = 0@f$, with @f$K=k^2@f$).
1224 
1225  The argument K_squared may be a real or a complex-valued vector.
1226 
1227  @ingroup asm
1228  */
1229  template<typename MAT, typename VECT>
1230  void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
1231  const mesh_fem &mf_data, const VECT &K_squared,
1232  const mesh_region &rg = mesh_region::all_convexes()) {
1233  typedef typename gmm::linalg_traits<MAT>::value_type T;
1234  asm_Helmholtz_(M, mim, mf_u, &mf_data, K_squared, rg, T());
1235  }
1236 
1237  template<typename MAT, typename VECT, typename T>
1238  void asm_Helmholtz_(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
1239  const mesh_fem *mf_data, const VECT &K_squared,
1240  const mesh_region &rg, T) {
1241  asm_real_or_complex_1_param_mat
1242  (M, mim, mf_u, mf_data, K_squared, rg,
1243  "(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u");
1244  }
1245 
1246  template<typename MAT, typename VECT, typename T>
1247  void asm_Helmholtz_(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
1248  const mesh_fem *mf_data, const VECT &K_squared,
1249  const mesh_region &rg, std::complex<T>) {
1250  // asm_real_or_complex_1_param_mat
1251  // (M, mim, mf_u, &mf_data, gmm::real_part(K_squared), rg,
1252  // "(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u");
1253  // asm_real_or_complex_1_param_mat
1254  // (M, mim, mf_u, &mf_data, gmm::imag_part(K_squared), rg,
1255  // "(A*Test_u).Test2_u");
1256 
1257 
1258  ga_workspace workspace;
1259  gmm::sub_interval Iur(0, mf_u.nb_dof()), Iui(Iur.last(), mf_u.nb_dof());
1260  base_vector u(mf_u.nb_dof());
1261  base_vector AR(gmm::vect_size(K_squared)), AI(gmm::vect_size(K_squared));
1262  gmm::copy(gmm::real_part(K_squared), AR);
1263  gmm::copy(gmm::imag_part(K_squared), AI);
1264  workspace.add_fem_variable("u", mf_u, Iur, u);
1265  workspace.add_fem_variable("ui", mf_u, Iui, u);
1266 
1267  if (mf_data) {
1268  workspace.add_fem_constant("A", *mf_data, AR);
1269  workspace.add_fem_constant("AI", *mf_data, AI);
1270  } else {
1271  workspace.add_fixed_size_constant("A", AR);
1272  workspace.add_fixed_size_constant("AI", AI);
1273  }
1274  workspace.add_expression("(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u",
1275  mim, rg);
1276  workspace.add_expression("(AI*Test_ui).Test2_ui", mim, rg);
1277  workspace.assembly(2);
1278 
1279  if (gmm::mat_nrows(workspace.assembled_matrix()))
1280  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(),Iur,Iur),
1281  const_cast<MAT &>(M));
1282  if (gmm::mat_nrows(workspace.assembled_matrix()) > mf_u.nb_dof())
1283  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(),Iui,Iui),
1284  gmm::imag_part(const_cast<MAT &>(M)));
1285 
1286 
1287 
1288  }
1289 
1290  /**
1291  assembly of the term @f$\int_\Omega Kuv - \nabla u.\nabla v@f$,
1292  for the helmholtz equation (@f$\Delta u + k^2u = 0@f$, with @f$K=k^2@f$).
1293 
1294  The argument K_squared may be a real or a complex-valued scalar.
1295 
1296  @ingroup asm
1297  */
1298  template<typename MAT, typename VECT>
1299  void asm_homogeneous_Helmholtz
1300  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const VECT &K_squared,
1301  const mesh_region &rg = mesh_region::all_convexes()) {
1302  typedef typename gmm::linalg_traits<MAT>::value_type T;
1303  asm_Helmholtz_(M, mim, mf_u, 0, K_squared, rg, T());
1304  }
1305 
1306  enum { ASMDIR_BUILDH = 1, ASMDIR_BUILDR = 2, ASMDIR_SIMPLIFY = 4,
1307  ASMDIR_BUILDALL = 7 };
1308 
1309  /**
1310  Assembly of Dirichlet constraints @f$ u(x) = r(x) @f$ in a weak form
1311  @f[ \int_{\Gamma} u(x)v(x) = \int_{\Gamma} r(x)v(x) \forall v@f],
1312  where @f$ v @f$ is in
1313  the space of multipliers corresponding to mf_mult.
1314 
1315  size(r_data) = Q * nb_dof(mf_rh);
1316 
1317  A simplification can be done when the fem for u and r are the same and
1318  when the fem for the multipliers is of same dimension as the one for u.
1319  version = |ASMDIR_BUILDH : build H
1320  |ASMDIR_BUILDR : build R
1321  |ASMDIR_SIMPLIFY : simplify
1322  |ASMDIR_BUILDALL : do everything.
1323 
1324  @ingroup asm
1325  */
1326 
1327  template<typename MAT, typename VECT1, typename VECT2>
1328  void asm_dirichlet_constraints
1329  (MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u,
1330  const mesh_fem &mf_mult, const mesh_fem &mf_r,
1331  const VECT2 &r_data, const mesh_region &region,
1332  int version = ASMDIR_BUILDALL) {
1333  typedef typename gmm::linalg_traits<VECT1>::value_type value_type;
1334  typedef typename gmm::number_traits<value_type>::magnitude_type magn_type;
1335 
1336  if ((version & ASMDIR_SIMPLIFY) &&
1337  (mf_u.is_reduced() || mf_mult.is_reduced() || mf_r.is_reduced())) {
1338  GMM_WARNING1("Sorry, no simplification for reduced fems");
1339  version = (version & (ASMDIR_BUILDR | ASMDIR_BUILDH));
1340  }
1341 
1342  region.from_mesh(mim.linked_mesh()).error_if_not_faces();
1343  GMM_ASSERT1(mf_r.get_qdim() == 1,
1344  "invalid data mesh fem (Qdim=1 required)");
1345  if (version & ASMDIR_BUILDH) {
1346  asm_mass_matrix(H, mim, mf_mult, mf_u, region);
1347  gmm::clean(H, gmm::default_tol(magn_type())
1348  * gmm::mat_maxnorm(H) * magn_type(1000));
1349  }
1350  if (version & ASMDIR_BUILDR)
1351  asm_source_term(R, mim, mf_mult, mf_r, r_data, region);
1352 
1353  // Verifications and simplifications
1354 
1355  pfem pf_u, pf_r, pf_m;
1356  bool warning_msg1 = false, warning_msg2 = false;
1357  dal::bit_vector simplifiable_dofs, nonsimplifiable_dofs;
1358  std::vector<size_type> simplifiable_indices(mf_mult.nb_basic_dof());
1359  std::vector<value_type> simplifiable_values(mf_mult.nb_basic_dof());
1360  std::vector<value_type> v1, v2, v3;
1361 
1362  for (mr_visitor v(region); !v.finished(); v.next()) {
1363  GMM_ASSERT1(v.is_face(), "attempt to impose a dirichlet "
1364  "on the interior of the domain!");
1365  size_type cv = v.cv(); short_type f = v.f();
1366 
1367  GMM_ASSERT1(mf_u.convex_index().is_in(cv) &&
1368  mf_r.convex_index().is_in(cv) &&
1369  mf_mult.convex_index().is_in(cv),
1370  "attempt to impose a dirichlet "
1371  "condition on a convex with no FEM!");
1372  pf_u = mf_u.fem_of_element(cv);
1373  pf_r = mf_r.fem_of_element(cv);
1374  pf_m = mf_mult.fem_of_element(cv);
1375 
1376  if (!pf_m->is_lagrange() && !warning_msg1) {
1377  GMM_WARNING3("Dirichlet condition with non-lagrange multiplier fem. "
1378  "see the documentation about Dirichlet conditions.");
1379  warning_msg1 = true;
1380  }
1381 
1382  if (!(version & ASMDIR_SIMPLIFY)) continue;
1383 
1384  mesh_fem::ind_dof_face_ct pf_u_ct
1385  = mf_u.ind_basic_dof_of_face_of_element(cv, f);
1386  mesh_fem::ind_dof_face_ct pf_r_ct
1387  = mf_r.ind_basic_dof_of_face_of_element(cv, f);
1388  mesh_fem::ind_dof_face_ct pf_m_ct
1389  = mf_mult.ind_basic_dof_of_face_of_element(cv, f);
1390 
1391  size_type pf_u_nbdf = pf_u_ct.size();
1392  size_type pf_m_nbdf = pf_m_ct.size();
1393  size_type pf_u_nbdf_loc = pf_u->structure(cv)->nb_points_of_face(f);
1394  size_type pf_m_nbdf_loc = pf_m->structure(cv)->nb_points_of_face(f);
1395  // size_type pf_r_nbdf_loc = pf_r->structure(cv)->nb_points_of_face(f);
1396 
1397  if (pf_u_nbdf < pf_m_nbdf && !warning_msg2) {
1398  GMM_WARNING2("Dirichlet condition with a too rich multiplier fem. "
1399  "see the documentation about Dirichlet conditions.");
1400  warning_msg2 = true;
1401  }
1402 
1403  if (pf_u != pf_r || pf_u_nbdf != pf_m_nbdf ||
1404  ((pf_u != pf_r) && (pf_u_nbdf_loc != pf_m_nbdf_loc))) {
1405  for (size_type i = 0; i < pf_m_nbdf; ++i)
1406  nonsimplifiable_dofs.add(pf_m_ct[i]);
1407  continue;
1408  }
1409 
1410  for (size_type i = 0; i < pf_m_nbdf; ++i) {
1411  simplifiable_dofs.add(pf_m_ct[i]);
1412  simplifiable_indices[pf_m_ct[i]] = pf_u_ct[i];
1413  }
1414 
1415  if (!(version & ASMDIR_BUILDR)) continue;
1416 
1417  if (pf_u == pf_r) { // simplest simplification.
1418  size_type Qratio = mf_u.get_qdim() / mf_r.get_qdim();
1419  for (size_type i = 0; i < pf_m_nbdf; ++i) {
1420  simplifiable_values[pf_m_ct[i]]
1421  = r_data[pf_r_ct[i/Qratio]*Qratio+(i%Qratio)];
1422  }
1423  }
1424  }
1425 
1426  if (version & ASMDIR_SIMPLIFY) {
1427  if (simplifiable_dofs.card() > 0)
1428  { GMM_TRACE3("Simplification of the Dirichlet condition"); }
1429  else
1430  GMM_TRACE3("Sorry, no simplification of the Dirichlet condition");
1431  if (nonsimplifiable_dofs.card() > 0 && simplifiable_dofs.card() > 0)
1432  GMM_WARNING3("Partial simplification of the Dirichlet condition");
1433 
1434  for (dal::bv_visitor i(simplifiable_dofs); !i.finished(); ++i)
1435  if (!(nonsimplifiable_dofs[i])) {
1436  if (version & ASMDIR_BUILDH) { /* "erase" the row i */
1437  const mesh::ind_cv_ct &cv_ct = mf_mult.convex_to_basic_dof(i);
1438  for (size_type j = 0; j < cv_ct.size(); ++j) {
1439  size_type cv = cv_ct[j];
1440  for (size_type k=0; k < mf_u.nb_basic_dof_of_element(cv); ++k)
1441  H(i, mf_u.ind_basic_dof_of_element(cv)[k]) = value_type(0);
1442  }
1443  H(i, simplifiable_indices[i]) = value_type(1);
1444  }
1445  if (version & ASMDIR_BUILDR) R[i] = simplifiable_values[i];
1446  }
1447  }
1448  }
1449 
1450  template<typename MATRM, typename VECT1, typename VECT2>
1451  void assembling_Dirichlet_condition
1452  (MATRM &RM, VECT1 &B, const getfem::mesh_fem &mf, size_type boundary,
1453  const VECT2 &F) {
1454  // Works only for Lagrange dofs.
1455  size_type Q=mf.get_qdim();
1456  GMM_ASSERT1(!(mf.is_reduced()), "This function is not adapted to "
1457  "reduced finite element methods");
1458  dal::bit_vector nndof = mf.basic_dof_on_region(boundary);
1459  getfem::pfem pf1;
1460  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
1461  pf1 = mf.fem_of_element(cv);
1462  getfem::pdof_description ldof = getfem::lagrange_dof(pf1->dim());
1463  size_type nbd = pf1->nb_dof(cv);
1464  for (size_type i = 0; i < nbd; i++) {
1465  size_type dof1 = mf.ind_basic_dof_of_element(cv)[i*Q];
1466  if (nndof.is_in(dof1) && pf1->dof_types()[i] == ldof) {
1467  // cout << "dof : " << i << endl;
1468  for (size_type j = 0; j < nbd; j++) {
1469  size_type dof2 = mf.ind_basic_dof_of_element(cv)[j*Q];
1470  for (size_type k = 0; k < Q; ++k)
1471  for (size_type l = 0; l < Q; ++l) {
1472  if (!(nndof.is_in(dof2)) &&
1473  getfem::dof_compatibility(pf1->dof_types()[j],
1474  getfem::lagrange_dof(pf1->dim())))
1475  B[dof2+k] -= RM(dof2+k, dof1+l) * F[dof1+l];
1476  RM(dof2+k, dof1+l) = RM(dof1+l, dof2+k) = 0;
1477  }
1478  }
1479  for (size_type k = 0; k < Q; ++k)
1480  { RM(dof1+k, dof1+k) = 1; B[dof1+k] = F[dof1+k]; }
1481  }
1482  }
1483  }
1484  }
1485 
1486  /**
1487  Assembly of generalized Dirichlet constraints h(x)u(x) = r(x),
1488  where h is a QxQ matrix field (Q == mf_u.get_qdim()), outputs a
1489  (under-determined) linear system HU=R.
1490 
1491  size(h_data) = Q^2 * nb_dof(mf_rh);
1492  size(r_data) = Q * nb_dof(mf_rh);
1493 
1494  This function tries hard to make H diagonal or mostly diagonal:
1495  this function is able to "simplify" the dirichlet constraints (see below)
1496  version = |ASMDIR_BUILDH : build H
1497  |ASMDIR_BUILDR : build R
1498  |ASMDIR_SIMPLIFY : simplify
1499  |ASMDIR_BUILDALL : do everything.
1500 
1501  @ingroup asm
1502  */
1503 
1504  template<typename MAT, typename VECT1, typename VECT2, typename VECT3>
1505  void asm_generalized_dirichlet_constraints
1506  (MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u,
1507  const mesh_fem &mf_h, const mesh_fem &mf_r, const VECT2 &h_data,
1508  const VECT3 &r_data, const mesh_region &region,
1509  int version = ASMDIR_BUILDALL) {
1510  pfem pf_u, pf_rh;
1511 
1512  if ((version & ASMDIR_SIMPLIFY) &&
1513  (mf_u.is_reduced() || mf_h.is_reduced() || mf_r.is_reduced())) {
1514  GMM_WARNING1("Sorry, no simplification for reduced fems");
1515  version = (version & ASMDIR_BUILDR);
1516  }
1517 
1518  region.from_mesh(mim.linked_mesh()).error_if_not_faces();
1519  GMM_ASSERT1(mf_h.get_qdim() == 1 && mf_r.get_qdim() == 1,
1520  "invalid data mesh fem (Qdim=1 required)");
1521  if (version & ASMDIR_BUILDH) {
1522  asm_qu_term(H, mim, mf_u, mf_h, h_data, region);
1523  std::vector<size_type> ind(0);
1524  dal::bit_vector bdof = mf_u.basic_dof_on_region(region);
1525  // gmm::clean(H, 1E-15 * gmm::mat_maxnorm(H));
1526  for (size_type i = 0; i < mf_u.nb_dof(); ++i)
1527  if (!(bdof[i])) ind.push_back(i);
1528  gmm::clear(gmm::sub_matrix(H, gmm::sub_index(ind)));
1529  }
1530  if (version & ASMDIR_BUILDR)
1531  asm_source_term(R, mim, mf_u, mf_r, r_data, region);
1532  if (!(version & ASMDIR_SIMPLIFY)) return;
1533 
1534  /* step 2 : simplification of simple dirichlet conditions */
1535  if (&mf_r == &mf_h) {
1536  for (mr_visitor v(region); !v.finished(); v.next()) {
1537  size_type cv = v.cv();
1538  short_type f = v.f();
1539 
1540  GMM_ASSERT1(mf_u.convex_index().is_in(cv) &&
1541  mf_r.convex_index().is_in(cv),
1542  "attempt to impose a dirichlet "
1543  "condition on a convex with no FEM!");
1544 
1545  if (f >= mf_u.linked_mesh().structure_of_convex(cv)->nb_faces())
1546  continue;
1547  pf_u = mf_u.fem_of_element(cv);
1548  pf_rh = mf_r.fem_of_element(cv);
1549  /* don't try anything with vector elements */
1550  if (mf_u.fem_of_element(cv)->target_dim() != 1) continue;
1551  bgeot::pconvex_structure cvs_u = pf_u->structure(cv);
1552  bgeot::pconvex_structure cvs_rh = pf_rh->structure(cv);
1553  for (size_type i = 0; i < cvs_u->nb_points_of_face(f); ++i) {
1554 
1555  size_type Q = mf_u.get_qdim(); // pf_u->target_dim() (==1)
1556 
1557  size_type ind_u = cvs_u->ind_points_of_face(f)[i];
1558  pdof_description tdof_u = pf_u->dof_types()[ind_u];
1559 
1560  for (size_type j = 0; j < cvs_rh->nb_points_of_face(f); ++j) {
1561  size_type ind_rh = cvs_rh->ind_points_of_face(f)[j];
1562  pdof_description tdof_rh = pf_rh->dof_types()[ind_rh];
1563  /*
1564  same kind of dof and same location of dof ?
1565  => then the previous was not useful for this dofs (introducing
1566  a mass matrix which is not diagonal in the constraints matrix)
1567  -> the constraint is simplified:
1568  we replace \int{(H_j.psi_j)*phi_i}=\int{R_j.psi_j} (sum over j)
1569  with H_j*phi_i = R_j
1570  --> Le principe peut être faux : non identique à la projection
1571  L^2 et peut entrer en conccurence avec les autres ddl -> a revoir
1572  */
1573  if (tdof_u == tdof_rh &&
1574  gmm::vect_dist2_sqr((*(pf_u->node_tab(cv)))[ind_u],
1575  (*(pf_rh->node_tab(cv)))[ind_rh])
1576  < 1.0E-14) {
1577  /* the dof might be "duplicated" */
1578  for (size_type q = 0; q < Q; ++q) {
1579  size_type dof_u = mf_u.ind_basic_dof_of_element(cv)[ind_u*Q+q];
1580 
1581  /* "erase" the row */
1582  if (version & ASMDIR_BUILDH)
1583  for (size_type k=0; k < mf_u.nb_basic_dof_of_element(cv); ++k)
1584  H(dof_u, mf_u.ind_basic_dof_of_element(cv)[k]) = 0.0;
1585 
1586  size_type dof_rh = mf_r.ind_basic_dof_of_element(cv)[ind_rh];
1587  /* set the "simplified" row */
1588  if (version & ASMDIR_BUILDH)
1589  for (unsigned jj=0; jj < Q; jj++) {
1590  size_type dof_u2
1591  = mf_u.ind_basic_dof_of_element(cv)[ind_u*Q+jj];
1592  H(dof_u, dof_u2) = h_data[(jj*Q+q) + Q*Q*(dof_rh)];
1593  }
1594  if (version & ASMDIR_BUILDR) R[dof_u] = r_data[dof_rh*Q+q];
1595  }
1596  }
1597  }
1598  }
1599  }
1600  }
1601  }
1602 
1603  /**
1604  Build an orthogonal basis of the kernel of H in NS, gives the
1605  solution of minimal norm of H*U = R in U0 and return the
1606  dimension of the kernel. The function is based on a
1607  Gramm-Schmidt algorithm.
1608 
1609  @ingroup asm
1610  */
1611  template<typename MAT1, typename MAT2, typename VECT1, typename VECT2>
1612  size_type Dirichlet_nullspace(const MAT1 &H, MAT2 &NS,
1613  const VECT1 &R, VECT2 &U0) {
1614 
1615  typedef typename gmm::linalg_traits<MAT1>::value_type T;
1616  typedef typename gmm::number_traits<T>::magnitude_type MAGT;
1617  typedef typename gmm::temporary_vector<MAT1>::vector_type TEMP_VECT;
1618  MAGT tol = gmm::default_tol(MAGT());
1619  MAGT norminfH = gmm::mat_maxnorm(H);
1620  size_type nbd = gmm::mat_ncols(H), nbase = 0, nbr = gmm::mat_nrows(H);
1621  TEMP_VECT aux(nbr), e(nbd), f(nbd);
1622  dal::dynamic_array<TEMP_VECT> base_img;
1623  dal::dynamic_array<TEMP_VECT> base_img_inv;
1624  size_type nb_bimg = 0;
1625  gmm::clear(NS);
1626 
1627  if (!(gmm::is_col_matrix(H)))
1628  GMM_WARNING2("Dirichlet_nullspace inefficient for a row matrix H");
1629  // First, detection of null columns of H, and already orthogonals
1630  // vectors of the image of H.
1631  dal::bit_vector nn;
1632  for (size_type i = 0; i < nbd; ++i) {
1633  gmm::clear(e); e[i] = T(1);
1634  gmm::mult(H, e, aux);
1635  MAGT n = gmm::vect_norm2(aux);
1636 
1637  if (n < norminfH*tol*1000) {
1638  NS(i, nbase++) = T(1); nn[i] = true;
1639  }
1640  else {
1641  bool good = true;
1642  for (size_type j = 0; j < nb_bimg; ++j)
1643  if (gmm::abs(gmm::vect_sp(aux, base_img[j])) > MAGT(0))
1644  { good = false; break; }
1645  if (good) {
1646  gmm::copy(e, f);
1647  gmm::scale(f, T(MAGT(1)/n)); gmm::scale(aux, T(MAGT(1)/n));
1648  base_img_inv[nb_bimg] = TEMP_VECT(nbd);
1649  gmm::copy(f, base_img_inv[nb_bimg]);
1650  gmm::clean(aux, gmm::vect_norminf(aux)*tol);
1651  base_img[nb_bimg] = TEMP_VECT(nbr);
1652  gmm::copy(aux, base_img[nb_bimg++]);
1653  nn[i] = true;
1654  }
1655  }
1656  }
1657  size_type nb_triv_base = nbase;
1658 
1659  for (size_type i = 0; i < nbd; ++i) {
1660  if (!(nn[i])) {
1661  gmm::clear(e); e[i] = T(1); gmm::clear(f); f[i] = T(1);
1662  gmm::mult(H, e, aux);
1663  for (size_type j = 0; j < nb_bimg; ++j) {
1664  T c = gmm::vect_sp(aux, base_img[j]);
1665  // if (gmm::abs(c > 1.0E-6) { // à scaler sur l'ensemble de H ...
1666  if (c != T(0)) {
1667  gmm::add(gmm::scaled(base_img[j], -c), aux);
1668  gmm::add(gmm::scaled(base_img_inv[j], -c), f);
1669  }
1670  }
1671  if (gmm::vect_norm2(aux) < norminfH*tol*MAGT(10000)) {
1672  gmm::copy(f, gmm::mat_col(NS, nbase++));
1673  }
1674  else {
1675  MAGT n = gmm::vect_norm2(aux);
1676  gmm::scale(f, T(MAGT(1)/n)); gmm::scale(aux, T(MAGT(1)/n));
1677  gmm::clean(f, tol*gmm::vect_norminf(f));
1678  gmm::clean(aux, tol*gmm::vect_norminf(aux));
1679  base_img_inv[nb_bimg] = TEMP_VECT(nbd);
1680  gmm::copy(f, base_img_inv[nb_bimg]);
1681  base_img[nb_bimg] = TEMP_VECT(nbr);
1682  gmm::copy(aux, base_img[nb_bimg++]);
1683  }
1684  }
1685  }
1686  // Compute a solution in U0
1687  gmm::clear(U0);
1688  for (size_type i = 0; i < nb_bimg; ++i) {
1689  T c = gmm::vect_sp(base_img[i], R);
1690  gmm::add(gmm::scaled(base_img_inv[i], c), U0);
1691  }
1692  // Orthogonalisation of the basis of the kernel of H.
1693  for (size_type i = nb_triv_base; i < nbase; ++i) {
1694  for (size_type j = nb_triv_base; j < i; ++j) {
1695  T c = gmm::vect_sp(gmm::mat_col(NS,i), gmm::mat_col(NS,j));
1696  if (c != T(0))
1697  gmm::add(gmm::scaled(gmm::mat_col(NS,j), -c), gmm::mat_col(NS,i));
1698  }
1699  gmm::scale(gmm::mat_col(NS,i),
1700  T(1) / gmm::vect_norm2(gmm::mat_col(NS,i)));
1701  }
1702  // projection of U0 on the orthogonal to the kernel.
1703  for (size_type j = nb_triv_base; j < nbase; ++j) {
1704  T c = gmm::vect_sp(gmm::mat_col(NS,j), U0);
1705  if (c != T(0))
1706  gmm::add(gmm::scaled(gmm::mat_col(NS,j), -c), U0);
1707  }
1708  // Test ...
1709  gmm::mult(H, U0, gmm::scaled(R, T(-1)), aux);
1710  if (gmm::vect_norm2(aux) > gmm::vect_norm2(U0)*tol*MAGT(10000))
1711  GMM_WARNING2("Dirichlet condition not well inverted: residu="
1712  << gmm::vect_norm2(aux));
1713 
1714  return nbase;
1715  }
1716 
1717 } /* end of namespace getfem. */
1718 
1719 
1720 #endif /* GETFEM_ASSEMBLING_H__ */
getfem::mesh_fem::nb_basic_dof_of_element
virtual size_type nb_basic_dof_of_element(size_type cv) const
Return the number of degrees of freedom attached to a given convex.
Definition: getfem_mesh_fem.h:493
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:560
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:66
getfem::asm_stiffness_matrix_for_homogeneous_scalar_elliptic_componentwise
void asm_stiffness_matrix_for_homogeneous_scalar_elliptic_componentwise(MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
The same but with a constant matrix.
Definition: getfem_assembling.h:1185
getfem::pdof_description
dof_description * pdof_description
Type representing a pointer on a dof_description.
Definition: getfem_fem.h:155
getfem::asm_stiffness_matrix_for_linear_elasticity
void asm_stiffness_matrix_for_linear_elasticity(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &LAMBDA, const VECT &MU, const mesh_region &rg=mesh_region::all_convexes())
Stiffness matrix for linear elasticity, with Lamé coefficients.
Definition: getfem_assembling.h:930
getfem::asm_stiffness_matrix_for_scalar_elliptic
void asm_stiffness_matrix_for_scalar_elliptic(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 , where is a (symmetric positive definite) NxN matrix.
Definition: getfem_assembling.h:1151
getfem::asm_stokes_B
void asm_stokes_B(const MAT &B, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_p, const mesh_region &rg=mesh_region::all_convexes())
Build the mixed pressure term .
Definition: getfem_assembling.h:1030
getfem::asm_H2_semi_norm
scalar_type asm_H2_semi_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute , U might be real or complex.
Definition: getfem_assembling.h:422
getfem::mesh_region::all_convexes
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
Definition: getfem_mesh_region.h:154
getfem::asm_H2_dist
scalar_type asm_H2_dist(const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1, const mesh_fem &mf2, const VEC2 &U2, const mesh_region &rg=mesh_region::all_convexes())
Compute the H2 distance between U1 and U2.
Definition: getfem_assembling.h:548
getfem::asm_homogeneous_Helmholtz
void asm_homogeneous_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const VECT &K_squared, const mesh_region &rg=mesh_region::all_convexes())
assembly of the term , for the helmholtz equation ( , with ).
Definition: getfem_assembling.h:1300
getfem::asm_stiffness_matrix_for_laplacian_componentwise
void asm_stiffness_matrix_for_laplacian_componentwise(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())
The same as getfem::asm_stiffness_matrix_for_laplacian , but on each component of mf when mf has a qd...
Definition: getfem_assembling.h:1121
getfem::asm_mass_matrix
void asm_mass_matrix(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_region &rg=mesh_region::all_convexes())
generic mass matrix assembly (on the whole mesh or on the specified convex set or boundary) ...
Definition: getfem_assembling.h:697
getfem::asm_stiffness_matrix_for_homogeneous_laplacian_componentwise
void asm_stiffness_matrix_for_homogeneous_laplacian_componentwise(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_region &rg=mesh_region::all_convexes())
assembly of .
Definition: getfem_assembling.h:1096
getfem::asm_stiffness_matrix_for_homogeneous_vector_elliptic
void asm_stiffness_matrix_for_homogeneous_vector_elliptic(MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
Assembly of , where is a NxNxQxQ (symmetric positive definite) constant tensor.
Definition: getfem_assembling.h:1213
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash! us...
Definition: getfem_mesh_fem.h:443
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:833
bgeot::mesh_structure::structure_of_convex
pconvex_structure structure_of_convex(size_type ic) const
Return the pconvex_structure of the convex ic.
Definition: bgeot_mesh_structure.h:112
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::asm_stiffness_matrix_for_homogeneous_scalar_elliptic
void asm_stiffness_matrix_for_homogeneous_scalar_elliptic(MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
The same but with a constant matrix.
Definition: getfem_assembling.h:1162
getfem::asm_stiffness_matrix_for_homogeneous_laplacian
void asm_stiffness_matrix_for_homogeneous_laplacian(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_region &rg=mesh_region::all_convexes())
assembly of .
Definition: getfem_assembling.h:1066
getfem::asm_dirichlet_constraints
void asm_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 &region, int version=ASMDIR_BUILDALL)
Assembly of Dirichlet constraints in a weak form , where is in the space of multipliers correspondi...
Definition: getfem_assembling.h:1329
getfem::asm_generalized_dirichlet_constraints
void asm_generalized_dirichlet_constraints(MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_h, const mesh_fem &mf_r, const VECT2 &h_data, const VECT3 &r_data, const mesh_region &region, int version=ASMDIR_BUILDALL)
Assembly of generalized Dirichlet constraints h(x)u(x) = r(x), where h is a QxQ matrix field (Q == mf...
Definition: getfem_assembling.h:1506
getfem::mesh_fem::convex_to_basic_dof
virtual const mesh::ind_cv_ct & convex_to_basic_dof(size_type d) const
Return the list of convexes attached to the specified dof.
Definition: getfem_mesh_fem.cc:267
getfem::asm_H1_dist
scalar_type asm_H1_dist(const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1, const mesh_fem &mf2, const VEC2 &U2, mesh_region rg=mesh_region::all_convexes())
Compute the H1 distance between U1 and U2.
Definition: getfem_assembling.h:355
getfem::mesh_region::from_mesh
const mesh_region & from_mesh(const mesh &m) const
For regions which have been built with just a number &#39;id&#39;, from_mesh(m) sets the current region to &#39;m...
Definition: getfem_mesh_region.cc:68
getfem::asm_stiffness_matrix_for_homogeneous_linear_elasticity
void asm_stiffness_matrix_for_homogeneous_linear_elasticity(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &LAMBDA, const VECT &MU, const mesh_region &rg=mesh_region::all_convexes())
Stiffness matrix for linear elasticity, with constant Lamé coefficients.
Definition: getfem_assembling.h:972
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
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:849
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::dof_compatibility
bool dof_compatibility(pdof_description, pdof_description)
Says if the two dofs can be identified.
Definition: getfem_fem.cc:597
getfem::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:555
getfem::mesh_fem::basic_dof_on_region
virtual dal::bit_vector basic_dof_on_region(const mesh_region &b) const
Get a list of basic dof lying on a given mesh_region.
Definition: getfem_mesh_fem.cc:77
getfem::asm_stiffness_matrix_for_scalar_elliptic_componentwise
void asm_stiffness_matrix_for_scalar_elliptic_componentwise(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())
The same but on each component of mf when mf has a qdim > 1.
Definition: getfem_assembling.h:1173
getfem::asm_H1_semi_dist
scalar_type asm_H1_semi_dist(const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1, const mesh_fem &mf2, const VEC2 &U2, mesh_region rg=mesh_region::all_convexes())
Compute the H1 semi-distance between U1 and U2, defined on two different mesh_fems (but sharing the s...
Definition: getfem_assembling.h:234
dal::dynamic_array
Dynamic Array.
Definition: dal_basic.h:47
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:311
getfem::asm_H2_norm
scalar_type asm_H2_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute the H2 norm of U (for C^1 elements).
Definition: getfem_assembling.h:534
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:252
getfem::Dirichlet_nullspace
size_type Dirichlet_nullspace(const MAT1 &H, MAT2 &NS, const VECT1 &R, VECT2 &U0)
Build an orthogonal basis of the kernel of H in NS, gives the solution of minimal norm of H*U = R in ...
Definition: getfem_assembling.h:1612
getfem_generic_assembly.h
A langage for generic assembly of pde boundary value problems.
getfem::mesh_fem::is_reduced
bool is_reduced() const
Return true if a reduction matrix is applied to the dofs.
Definition: getfem_mesh_fem.h:195
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_assembling.h:48
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:385
getfem::asm_stiffness_matrix_for_linear_elasticity_Hooke
void asm_stiffness_matrix_for_linear_elasticity_Hooke(MAT &RM, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &H, const mesh_region &rg=mesh_region::all_convexes())
Stiffness matrix for linear elasticity, with a general Hooke tensor.
Definition: getfem_assembling.h:1019
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:239
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:79
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:861
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:875
getfem::asm_L2_dist
scalar_type asm_L2_dist(const mesh_im &mim, const mesh_fem &mf1, const VEC1 &U1, const mesh_fem &mf2, const VEC2 &U2, mesh_region rg=mesh_region::all_convexes())
Compute the distance between U1 and U2, defined on two different mesh_fems (but sharing the same mesh...
Definition: getfem_assembling.h:116
getfem::asm_mass_matrix_homogeneous_param
void asm_mass_matrix_homogeneous_param(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const VECT &F, const mesh_region &rg=mesh_region::all_convexes())
generic mass matrix assembly with an additional constant parameter (on the whole mesh or on the speci...
Definition: getfem_assembling.h:822
getfem::asm_stiffness_matrix_for_vector_elliptic
void asm_stiffness_matrix_for_vector_elliptic(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 , where is a NxNxQxQ (symmetric positive definite) tensor defined on mf_data...
Definition: getfem_assembling.h:1199
getfem::asm_L2_norm
scalar_type asm_L2_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute , U might be real or complex
Definition: getfem_assembling.h:56
getfem::mesh_fem::ind_basic_dof_of_face_of_element
virtual ind_dof_face_ct ind_basic_dof_of_face_of_element(size_type cv, short_type f) const
Give an array of the dof numbers lying of a convex face (all degrees of freedom whose associated base...
Definition: getfem_mesh_fem.h:467
getfem::asm_stiffness_matrix_for_laplacian
void asm_stiffness_matrix_for_laplacian(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 , where is scalar.
Definition: getfem_assembling.h:1108
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:278
getfem::mesh_im::linked_mesh
const mesh & linked_mesh() const
Give a reference to the linked mesh of type mesh.
Definition: getfem_mesh_im.h:79
getfem_assembling_tensors.h
Generic assembly implementation.
getfem::mesh_fem::convex_index
const dal::bit_vector & convex_index() const
Get the set of convexes where a finite element has been assigned.
Definition: getfem_mesh_fem.h:191
getfem::asm_Helmholtz
void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_data, const VECT &K_squared, const mesh_region &rg=mesh_region::all_convexes())
assembly of the term , for the helmholtz equation ( , with ).
Definition: getfem_assembling.h:1230
getfem::asm_H1_norm
scalar_type asm_H1_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute the H1 norm of U.
Definition: getfem_assembling.h:302
getfem::mesh_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:449
getfem::asm_mass_matrix_param
void asm_mass_matrix_param(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_fem &mf2, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
generic mass matrix assembly with an additional parameter (on the whole mesh or on the specified boun...
Definition: getfem_assembling.h:765
getfem::asm_qu_term
void asm_qu_term(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_d, const VECT &Q, const mesh_region &rg)
assembly of
Definition: getfem_assembling.h:900
getfem::asm_H1_semi_norm
scalar_type asm_H1_semi_norm(const mesh_im &mim, const mesh_fem &mf, const VEC &U, const mesh_region &rg=mesh_region::all_convexes())
compute , U might be real or complex
Definition: getfem_assembling.h:179