GetFEM  5.4.3
getfem_assembling.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2000-2020 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  lumped mass matrix assembly from consistent mass matrix
830  */
831  template<typename MAT>
832  inline void asm_lumped_mass_matrix_for_first_order_from_consistent
833  (const MAT &M) {
834  size_type nbd = gmm::mat_ncols(M), nbr = gmm::mat_nrows(M);
835  GMM_ASSERT1(nbd == nbr, "mass matrix is not square");
836  typedef typename gmm::linalg_traits<MAT>::value_type T;
837  std::vector<T> V(nbd), W(nbr);
838  gmm::fill(V, T(1));
839  gmm::mult(M, V, W);
840  gmm::clear(const_cast<MAT &>(M));
841  for (size_type i =0; i < nbd; ++i) {
842  (const_cast<MAT &>(M))(i, i) = W[i];
843  }
844  }
845 
846  /**
847  lumped mass matrix assembly (on the whole mesh or on the specified
848  boundary)
849  @ingroup asm
850  */
851  template<typename MAT>
852  inline void asm_lumped_mass_matrix_for_first_order
853  (const MAT &M, const mesh_im &mim, const mesh_fem &mf1,
854  const mesh_region &rg = mesh_region::all_convexes()) {
855  asm_mass_matrix(M, mim, mf1, rg);
856  asm_lumped_mass_matrix_for_first_order_from_consistent(M);
857  }
858 
859  /**
860  lumped mass matrix assembly with an additional parameter
861  (on the whole mesh or on the specified boundary)
862  @ingroup asm
863  */
864  template<typename MAT, typename VECT>
865  inline void asm_lumped_mass_matrix_for_first_order_param
866  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_data,
867  const VECT &F, const mesh_region &rg = mesh_region::all_convexes()) {
868  asm_mass_matrix_param(M, mim, mf_u, mf_data, F, rg);
869  asm_lumped_mass_matrix_for_first_order_from_consistent(M);
870  }
871 
872  /**
873  source term (for both volumic sources and boundary (Neumann) sources).
874  @ingroup asm
875  */
876  template<typename VECT1, typename VECT2>
877  void asm_source_term(const VECT1 &B, const mesh_im &mim, const mesh_fem &mf,
878  const mesh_fem &mf_data, const VECT2 &F,
879  const mesh_region &rg = mesh_region::all_convexes()) {
880  GMM_ASSERT1(mf_data.get_qdim() == 1 ||
881  mf_data.get_qdim() == mf.get_qdim(),
882  "invalid data mesh fem (same Qdim or Qdim=1 required)");
883  asm_real_or_complex_1_param_vec
884  (const_cast<VECT1 &>(B), mim, mf, &mf_data, F, rg, "A:Test_u");
885  }
886 
887  /**
888  source term (for both volumic sources and boundary (Neumann) sources).
889  For an homogeneous term.
890  @ingroup asm
891  */
892  template<typename VECT1, typename VECT2>
893  void asm_homogeneous_source_term(const VECT1 &B, const mesh_im &mim,
894  const mesh_fem &mf, const VECT2 &F,
895  const mesh_region &rg = mesh_region::all_convexes()) {
896  asm_real_or_complex_1_param_vec
897  (const_cast<VECT1 &>(B), mim, mf, 0, F, rg, "A:Test_u");
898  }
899 
900  /**
901  Normal source term (for boundary (Neumann) condition).
902  @ingroup asm
903  */
904  template<typename VECT1, typename VECT2>
905  void asm_normal_source_term(VECT1 &B, const mesh_im &mim,
906  const mesh_fem &mf,
907  const mesh_fem &mf_data, const VECT2 &F,
908  const mesh_region &rg) {
909  asm_real_or_complex_1_param_vec(B, mim, mf, &mf_data, F, rg,
910  "(Reshape(A, qdim(u), meshdim).Normal):Test_u");
911  }
912 
913  /**
914  Homogeneous normal source term (for boundary (Neumann) condition).
915  @ingroup asm
916  */
917  template<typename VECT1, typename VECT2>
918  void asm_homogeneous_normal_source_term
919  (VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const VECT2 &F,
920  const mesh_region &rg)
921  { asm_real_or_complex_1_param_vec(B, mim, mf, 0,F,rg,
922  "(Reshape(A, qdim(u), meshdim).Normal):Test_u"); }
923 
924  /**
925  assembly of @f$\int{qu.v}@f$
926 
927  (if @f$u@f$ is a vector field of size @f$N@f$, @f$q@f$ is a square
928  matrix @f$N\times N@f$ used by assem_general_boundary_conditions
929 
930  convention: Q is of the form
931  Q1_11 Q2_11 ..... Qn_11
932  Q1_21 Q2_21 ..... Qn_21
933  Q1_12 Q2_12 ..... Qn_12
934  Q1_22 Q2_22 ..... Qn_22
935  if N = 2, and mf_d has n/N degree of freedom
936 
937  Q is a vector, so the matrix is assumed to be stored by columns
938  (fortran style)
939 
940  Works for both volumic assembly and boundary assembly
941  @ingroup asm
942  */
943  template<typename MAT, typename VECT>
944  void asm_qu_term(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
945  const mesh_fem &mf_d, const VECT &Q,
946  const mesh_region &rg) {
947  if (mf_d.get_qdim() == 1 && gmm::vect_size(Q) > mf_d.nb_dof())
948  asm_real_or_complex_1_param_mat
949  (M, mim,mf_u,&mf_d,Q,rg,"(Reshape(A,qdim(u),qdim(u)).Test_u):Test2_u");
950  else if (mf_d.get_qdim() == mf_u.get_qdim())
951  asm_real_or_complex_1_param_mat
952  (M, mim, mf_u, &mf_d, Q, rg, "(A*Test_u):Test2_u");
953  else GMM_ASSERT1(false, "invalid data mesh fem");
954  }
955 
956  template<typename MAT, typename VECT>
957  void asm_homogeneous_qu_term(MAT &M, const mesh_im &mim,
958  const mesh_fem &mf_u, const VECT &Q,
959  const mesh_region &rg) {
960  if (gmm::vect_size(Q) == 1)
961  asm_real_or_complex_1_param_mat
962  (M, mim,mf_u,0,Q,rg,"(A*Test_u):Test2_u");
963  else
964  asm_real_or_complex_1_param_mat
965  (M, mim,mf_u,0,Q,rg,"(Reshape(A,qdim(u),qdim(u)).Test_u):Test2_u");
966  }
967 
968  /**
969  Stiffness matrix for linear elasticity, with Lamé coefficients
970  @ingroup asm
971  */
972  template<class MAT, class VECT>
973  inline void asm_stiffness_matrix_for_linear_elasticity
974  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
975  const mesh_fem &mf_data, const VECT &LAMBDA, const VECT &MU,
976  const mesh_region &rg = mesh_region::all_convexes()) {
977  ga_workspace workspace;
978  gmm::sub_interval Iu(0, mf.nb_dof());
979  base_vector u(mf.nb_dof()), lambda(gmm::vect_size(LAMBDA));
980  base_vector mu(gmm::vect_size(MU));
981  gmm::copy(LAMBDA, lambda); gmm::copy(MU, mu);
982  workspace.add_fem_variable("u", mf, Iu, u);
983  workspace.add_fem_constant("lambda", mf_data, lambda);
984  workspace.add_fem_constant("mu", mf_data, mu);
985  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
986  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
987  workspace.assembly(2);
988  if (gmm::mat_nrows(workspace.assembled_matrix()))
989  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
990  }
991 
992  inline void asm_stiffness_matrix_for_linear_elasticity
993  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf,
994  const mesh_fem &mf_data, const model_real_plain_vector &LAMBDA,
995  const model_real_plain_vector &MU,
996  const mesh_region &rg = mesh_region::all_convexes()) {
997  ga_workspace workspace;
998  gmm::sub_interval Iu(0, mf.nb_dof());
999  base_vector u(mf.nb_dof());
1000  workspace.add_fem_variable("u", mf, Iu, u);
1001  workspace.add_fem_constant("lambda", mf_data, LAMBDA);
1002  workspace.add_fem_constant("mu", mf_data, MU);
1003  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
1004  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
1005  workspace.set_assembled_matrix(M);
1006  workspace.assembly(2);
1007  }
1008 
1009 
1010  /**
1011  Stiffness matrix for linear elasticity, with constant Lamé coefficients
1012  @ingroup asm
1013  */
1014  template<class MAT, class VECT>
1015  inline void asm_stiffness_matrix_for_homogeneous_linear_elasticity
1016  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
1017  const VECT &LAMBDA, const VECT &MU,
1018  const mesh_region &rg = mesh_region::all_convexes()) {
1019  ga_workspace workspace;
1020  gmm::sub_interval Iu(0, mf.nb_dof());
1021  base_vector u(mf.nb_dof()), lambda(gmm::vect_size(LAMBDA));
1022  base_vector mu(gmm::vect_size(MU));
1023  gmm::copy(LAMBDA, lambda); gmm::copy(MU, mu);
1024  workspace.add_fem_variable("u", mf, Iu, u);
1025  workspace.add_fixed_size_constant("lambda", lambda);
1026  workspace.add_fixed_size_constant("mu", mu);
1027  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
1028  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
1029  workspace.assembly(2);
1030  if (gmm::mat_nrows(workspace.assembled_matrix()))
1031  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
1032  }
1033 
1034 
1035  inline void asm_stiffness_matrix_for_homogeneous_linear_elasticity
1036  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf,
1037  const model_real_plain_vector &LAMBDA, const model_real_plain_vector &MU,
1038  const mesh_region &rg = mesh_region::all_convexes()) {
1039  ga_workspace workspace;
1040  gmm::sub_interval Iu(0, mf.nb_dof());
1041  base_vector u(mf.nb_dof());
1042  workspace.add_fem_variable("u", mf, Iu, u);
1043  workspace.add_fixed_size_constant("lambda", LAMBDA);
1044  workspace.add_fixed_size_constant("mu", MU);
1045  workspace.add_expression("((lambda*Div_Test_u)*Id(meshdim)"
1046  "+(2*mu)*Sym(Grad_Test_u)):Grad_Test2_u", mim, rg);
1047  workspace.set_assembled_matrix(M);
1048  workspace.assembly(2);
1049  }
1050 
1051  /**
1052  Stiffness matrix for linear elasticity, with a general Hooke
1053  tensor. This is more a demonstration of generic assembly than
1054  something useful !
1055 
1056  Note that this function is just an alias for
1057  asm_stiffness_matrix_for_vector_elliptic.
1058 
1059  @ingroup asm
1060  */
1061  template<typename MAT, typename VECT> void
1062  asm_stiffness_matrix_for_linear_elasticity_Hooke
1063  (MAT &RM, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1064  const VECT &H, const mesh_region &rg = mesh_region::all_convexes()) {
1065  asm_stiffness_matrix_for_vector_elliptic(RM, mim, mf, mf_data, H, rg);
1066  }
1067 
1068  /**
1069  Build the mixed pressure term @f$ B = - \int p.div u @f$
1070 
1071  @ingroup asm
1072  */
1073  template<typename MAT>
1074  inline void asm_stokes_B(const MAT &B, const mesh_im &mim,
1075  const mesh_fem &mf_u, const mesh_fem &mf_p,
1076  const mesh_region &rg=mesh_region::all_convexes()) {
1077  ga_workspace workspace;
1078  gmm::sub_interval Iu(0, mf_u.nb_dof()), Ip(Iu.last(), mf_p.nb_dof());
1079  base_vector u(mf_u.nb_dof()), p(mf_p.nb_dof());
1080  workspace.add_fem_variable("u", mf_u, Iu, u);
1081  workspace.add_fem_variable("p", mf_p, Ip, p);
1082  workspace.add_expression("Test_p*Div_Test2_u", mim, rg);
1083  workspace.assembly(2);
1084  if (gmm::mat_nrows(workspace.assembled_matrix()))
1085  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(), Ip, Iu),
1086  const_cast<MAT &>(B));
1087  }
1088 
1089  inline void asm_stokes_B(model_real_sparse_matrix &B, const mesh_im &mim,
1090  const mesh_fem &mf_u, const mesh_fem &mf_p,
1091  const mesh_region &rg=mesh_region::all_convexes()) {
1092  ga_workspace workspace;
1093  gmm::sub_interval Iu(0, mf_u.nb_dof()), Ip(0, mf_p.nb_dof());
1094  base_vector u(mf_u.nb_dof()), p(mf_p.nb_dof());
1095  workspace.add_fem_variable("u", mf_u, Iu, u);
1096  workspace.add_fem_variable("p", mf_p, Ip, p);
1097  workspace.add_expression("Test_p*Div_Test2_u", mim, rg);
1098  workspace.set_assembled_matrix(B);
1099  workspace.assembly(2);
1100  }
1101 
1102 
1103  /**
1104  assembly of @f$\int_\Omega \nabla u.\nabla v@f$.
1105 
1106  @ingroup asm
1107  */
1108  template<typename MAT>
1109  inline void asm_stiffness_matrix_for_homogeneous_laplacian
1110  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
1111  const mesh_region &rg = mesh_region::all_convexes()) {
1112  ga_workspace workspace;
1113  gmm::sub_interval Iu(0, mf.nb_dof());
1114  base_vector u(mf.nb_dof());
1115  workspace.add_fem_variable("u", mf, Iu, u);
1116  workspace.add_expression("Grad_Test_u:Grad_Test2_u", mim, rg);
1117  workspace.assembly(2);
1118  if (gmm::mat_nrows(workspace.assembled_matrix()))
1119  gmm::add(workspace.assembled_matrix(), const_cast<MAT &>(M));
1120  }
1121 
1122  inline void asm_stiffness_matrix_for_homogeneous_laplacian
1123  (model_real_sparse_matrix &M, const mesh_im &mim, const mesh_fem &mf,
1124  const mesh_region &rg = mesh_region::all_convexes()) {
1125  ga_workspace workspace;
1126  gmm::sub_interval Iu(0, mf.nb_dof());
1127  base_vector u(mf.nb_dof());
1128  workspace.add_fem_variable("u", mf, Iu, u);
1129  workspace.add_expression("Grad_Test_u:Grad_Test2_u", mim, rg);
1130  workspace.set_assembled_matrix(M);
1131  workspace.assembly(2);
1132  }
1133 
1134  /**
1135  assembly of @f$\int_\Omega \nabla u.\nabla v@f$.
1136  @ingroup asm
1137  */
1138  template<typename MAT>
1139  inline void asm_stiffness_matrix_for_homogeneous_laplacian_componentwise
1140  (const MAT &M, const mesh_im &mim, const mesh_fem &mf,
1141  const mesh_region &rg = mesh_region::all_convexes()) {
1142  return asm_stiffness_matrix_for_homogeneous_laplacian(M, mim, mf, rg);
1143  }
1144 
1145  /**
1146  assembly of @f$\int_\Omega a(x)\nabla u.\nabla v@f$ , where @f$a(x)@f$
1147  is scalar.
1148  @ingroup asm
1149  */
1150  template<typename MAT, typename VECT>
1151  inline void asm_stiffness_matrix_for_laplacian
1152  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1153  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1154  GMM_ASSERT1(mf_data.get_qdim() == 1
1155  && gmm::vect_size(A) == mf_data.nb_dof(), "invalid data");
1156  asm_real_or_complex_1_param_mat
1157  (M, mim, mf, &mf_data, A, rg, "(A*Grad_Test_u):Grad_Test2_u");
1158  }
1159 
1160  /** The same as getfem::asm_stiffness_matrix_for_laplacian , but on
1161  each component of mf when mf has a qdim > 1
1162  */
1163  template<typename MAT, typename VECT>
1164  inline void asm_stiffness_matrix_for_laplacian_componentwise
1165  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1166  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1167  asm_stiffness_matrix_for_laplacian(M, mim, mf, mf_data, A, rg);
1168  }
1169 
1170  /**
1171  assembly of @f$\int_\Omega A(x)\nabla u.\nabla v@f$, where @f$A(x)@f$
1172  is a (symmetric positive definite) NxN matrix.
1173  Arguments:
1174  @param M a sparse matrix of dimensions mf.nb_dof() x mf.nb_dof()
1175 
1176  @param mim the mesh_im.
1177 
1178  @param mf : the mesh_fem that describes the solution, with
1179  @c mf.get_qdim() == @c N.
1180 
1181  @param mf_data the mesh_fem that describes the coefficients of @c A
1182  (@c mf_data.get_qdim() == 1).
1183 
1184  @param A a (very large) vector, which is a flattened (n x n x
1185  mf_data.nb_dof()) 3D array. For each dof of mf_data, it contains
1186  the n x n coefficients of @f$A@f$. As usual, the order is the
1187  "fortran-order", i.e. @c A = [A_11(dof1) A_21(dof1) A_31(dof1)
1188  A_12(dof1) A_22(dof1) ... A_33(dof) A_11(dof2)
1189  .... A_33(lastdof)]
1190 
1191  @ingroup asm
1192  */
1193  template<typename MAT, typename VECT>
1194  void asm_stiffness_matrix_for_scalar_elliptic
1195  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1196  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1197  asm_real_or_complex_1_param_mat
1198  (M, mim, mf, &mf_data, A, rg,
1199  "(Reshape(A,meshdim,meshdim)*Grad_Test_u):Grad_Test2_u");
1200  }
1201 
1202  /** The same but with a constant matrix
1203  */
1204  template<typename MAT, typename VECT>
1205  void asm_stiffness_matrix_for_homogeneous_scalar_elliptic
1206  (MAT &M, const mesh_im &mim, const mesh_fem &mf,
1207  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1208  asm_real_or_complex_1_param_mat
1209  (M, mim, mf, 0, A, rg,
1210  "(Reshape(A,meshdim,meshdim)*Grad_Test_u):Grad_Test2_u");
1211  }
1212 
1213  /** The same but on each component of mf when mf has a qdim > 1
1214  */
1215  template<typename MAT, typename VECT>
1216  void asm_stiffness_matrix_for_scalar_elliptic_componentwise
1217  (MAT &M, const mesh_im &mim, const mesh_fem &mf,
1218  const mesh_fem &mf_data, const VECT &A,
1219  const mesh_region &rg = mesh_region::all_convexes()) {
1220  asm_real_or_complex_1_param_mat
1221  (M, mim, mf, &mf_data, A, rg,
1222  "(Grad_Test_u*(Reshape(A,meshdim,meshdim)')):Grad_Test2_u");
1223  }
1224 
1225  /** The same but with a constant matrix
1226  */
1227  template<typename MAT, typename VECT>
1228  void asm_stiffness_matrix_for_homogeneous_scalar_elliptic_componentwise
1229  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const VECT &A,
1230  const mesh_region &rg = mesh_region::all_convexes()) {
1231  asm_real_or_complex_1_param_mat
1232  (M, mim, mf, 0, A, rg,
1233  "(Grad_Test_u*(Reshape(A,meshdim,meshdim)')):Grad_Test2_u");
1234  }
1235 
1236 
1237  /**
1238  Assembly of @f$\int_\Omega A(x)\nabla u.\nabla v@f$, where @f$A(x)@f$
1239  is a NxNxQxQ (symmetric positive definite) tensor defined on mf_data.
1240  */
1241  template<typename MAT, typename VECT> void
1242  asm_stiffness_matrix_for_vector_elliptic
1243  (MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data,
1244  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1245  /* M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v) */
1246  asm_real_or_complex_1_param_mat
1247  (M, mim, mf, &mf_data, A, rg,
1248  "(Reshape(A,qdim(u),meshdim,qdim(u),meshdim):Grad_Test_u):Grad_Test2_u");
1249  }
1250 
1251  /**
1252  Assembly of @f$\int_\Omega A(x)\nabla u.\nabla v@f$, where @f$A(x)@f$
1253  is a NxNxQxQ (symmetric positive definite) constant tensor.
1254  */
1255  template<typename MAT, typename VECT> void
1256  asm_stiffness_matrix_for_homogeneous_vector_elliptic
1257  (MAT &M, const mesh_im &mim, const mesh_fem &mf,
1258  const VECT &A, const mesh_region &rg = mesh_region::all_convexes()) {
1259  /* M = a_{i,j,k,l}D_{i,j}(u)D_{k,l}(v) */
1260  asm_real_or_complex_1_param_mat
1261  (M, mim, mf, 0, A, rg,
1262  "(Reshape(A,qdim(u),meshdim,qdim(u),meshdim):Grad_Test_u):Grad_Test2_u");
1263  }
1264 
1265  /**
1266  assembly of the term @f$\int_\Omega Kuv - \nabla u.\nabla v@f$,
1267  for the helmholtz equation (@f$\Delta u + k^2u = 0@f$, with @f$K=k^2@f$).
1268 
1269  The argument K_squared may be a real or a complex-valued vector.
1270 
1271  @ingroup asm
1272  */
1273  template<typename MAT, typename VECT>
1274  void asm_Helmholtz(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
1275  const mesh_fem &mf_data, const VECT &K_squared,
1276  const mesh_region &rg = mesh_region::all_convexes()) {
1277  typedef typename gmm::linalg_traits<MAT>::value_type T;
1278  asm_Helmholtz_(M, mim, mf_u, &mf_data, K_squared, rg, T());
1279  }
1280 
1281  template<typename MAT, typename VECT, typename T>
1282  void asm_Helmholtz_(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
1283  const mesh_fem *mf_data, const VECT &K_squared,
1284  const mesh_region &rg, T) {
1285  asm_real_or_complex_1_param_mat
1286  (M, mim, mf_u, mf_data, K_squared, rg,
1287  "(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u");
1288  }
1289 
1290  template<typename MAT, typename VECT, typename T>
1291  void asm_Helmholtz_(MAT &M, const mesh_im &mim, const mesh_fem &mf_u,
1292  const mesh_fem *mf_data, const VECT &K_squared,
1293  const mesh_region &rg, std::complex<T>) {
1294  // asm_real_or_complex_1_param_mat
1295  // (M, mim, mf_u, &mf_data, gmm::real_part(K_squared), rg,
1296  // "(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u");
1297  // asm_real_or_complex_1_param_mat
1298  // (M, mim, mf_u, &mf_data, gmm::imag_part(K_squared), rg,
1299  // "(A*Test_u).Test2_u");
1300 
1301 
1302  ga_workspace workspace;
1303  gmm::sub_interval Iur(0, mf_u.nb_dof()), Iui(Iur.last(), mf_u.nb_dof());
1304  base_vector u(mf_u.nb_dof());
1305  base_vector AR(gmm::vect_size(K_squared)), AI(gmm::vect_size(K_squared));
1306  gmm::copy(gmm::real_part(K_squared), AR);
1307  gmm::copy(gmm::imag_part(K_squared), AI);
1308  workspace.add_fem_variable("u", mf_u, Iur, u);
1309  workspace.add_fem_variable("ui", mf_u, Iui, u);
1310 
1311  if (mf_data) {
1312  workspace.add_fem_constant("A", *mf_data, AR);
1313  workspace.add_fem_constant("AI", *mf_data, AI);
1314  } else {
1315  workspace.add_fixed_size_constant("A", AR);
1316  workspace.add_fixed_size_constant("AI", AI);
1317  }
1318  workspace.add_expression("(A*Test_u).Test2_u - Grad_Test_u:Grad_Test2_u",
1319  mim, rg);
1320  workspace.add_expression("(AI*Test_ui).Test2_ui", mim, rg);
1321  workspace.assembly(2);
1322 
1323  if (gmm::mat_nrows(workspace.assembled_matrix()))
1324  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(),Iur,Iur),
1325  const_cast<MAT &>(M));
1326  if (gmm::mat_nrows(workspace.assembled_matrix()) > mf_u.nb_dof())
1327  gmm::add(gmm::sub_matrix(workspace.assembled_matrix(),Iui,Iui),
1328  gmm::imag_part(const_cast<MAT &>(M)));
1329 
1330 
1331 
1332  }
1333 
1334  /**
1335  assembly of the term @f$\int_\Omega Kuv - \nabla u.\nabla v@f$,
1336  for the helmholtz equation (@f$\Delta u + k^2u = 0@f$, with @f$K=k^2@f$).
1337 
1338  The argument K_squared may be a real or a complex-valued scalar.
1339 
1340  @ingroup asm
1341  */
1342  template<typename MAT, typename VECT>
1343  void asm_homogeneous_Helmholtz
1344  (MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const VECT &K_squared,
1345  const mesh_region &rg = mesh_region::all_convexes()) {
1346  typedef typename gmm::linalg_traits<MAT>::value_type T;
1347  asm_Helmholtz_(M, mim, mf_u, 0, K_squared, rg, T());
1348  }
1349 
1350  enum { ASMDIR_BUILDH = 1, ASMDIR_BUILDR = 2, ASMDIR_SIMPLIFY = 4,
1351  ASMDIR_BUILDALL = 7 };
1352 
1353  /**
1354  Assembly of Dirichlet constraints @f$ u(x) = r(x) @f$ in a weak form
1355  @f[ \int_{\Gamma} u(x)v(x) = \int_{\Gamma} r(x)v(x) \forall v@f],
1356  where @f$ v @f$ is in
1357  the space of multipliers corresponding to mf_mult.
1358 
1359  size(r_data) = Q * nb_dof(mf_rh);
1360 
1361  A simplification can be done when the fem for u and r are the same and
1362  when the fem for the multipliers is of same dimension as the one for u.
1363  version = |ASMDIR_BUILDH : build H
1364  |ASMDIR_BUILDR : build R
1365  |ASMDIR_SIMPLIFY : simplify
1366  |ASMDIR_BUILDALL : do everything.
1367 
1368  @ingroup asm
1369  */
1370 
1371  template<typename MAT, typename VECT1, typename VECT2>
1372  void asm_dirichlet_constraints
1373  (MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u,
1374  const mesh_fem &mf_mult, const mesh_fem &mf_r,
1375  const VECT2 &r_data, const mesh_region &region,
1376  int version = ASMDIR_BUILDALL) {
1377  typedef typename gmm::linalg_traits<VECT1>::value_type value_type;
1378  typedef typename gmm::number_traits<value_type>::magnitude_type magn_type;
1379 
1380  if ((version & ASMDIR_SIMPLIFY) &&
1381  (mf_u.is_reduced() || mf_mult.is_reduced() || mf_r.is_reduced())) {
1382  GMM_WARNING1("Sorry, no simplification for reduced fems");
1383  version = (version & (ASMDIR_BUILDR | ASMDIR_BUILDH));
1384  }
1385 
1386  region.from_mesh(mim.linked_mesh()).error_if_not_faces();
1387  GMM_ASSERT1(mf_r.get_qdim() == 1,
1388  "invalid data mesh fem (Qdim=1 required)");
1389  if (version & ASMDIR_BUILDH) {
1390  asm_mass_matrix(H, mim, mf_mult, mf_u, region);
1391  gmm::clean(H, gmm::default_tol(magn_type())
1392  * gmm::mat_maxnorm(H) * magn_type(1000));
1393  }
1394  if (version & ASMDIR_BUILDR)
1395  asm_source_term(R, mim, mf_mult, mf_r, r_data, region);
1396 
1397  // Verifications and simplifications
1398 
1399  pfem pf_u, pf_r, pf_m;
1400  bool warning_msg1 = false, warning_msg2 = false;
1401  dal::bit_vector simplifiable_dofs, nonsimplifiable_dofs;
1402  std::vector<size_type> simplifiable_indices(mf_mult.nb_basic_dof());
1403  std::vector<value_type> simplifiable_values(mf_mult.nb_basic_dof());
1404  std::vector<value_type> v1, v2, v3;
1405 
1406  for (mr_visitor v(region); !v.finished(); v.next()) {
1407  GMM_ASSERT1(v.is_face(), "attempt to impose a dirichlet "
1408  "on the interior of the domain!");
1409  size_type cv = v.cv(); short_type f = v.f();
1410 
1411  GMM_ASSERT1(mf_u.convex_index().is_in(cv) &&
1412  mf_r.convex_index().is_in(cv) &&
1413  mf_mult.convex_index().is_in(cv),
1414  "attempt to impose a dirichlet "
1415  "condition on a convex with no FEM!");
1416  pf_u = mf_u.fem_of_element(cv);
1417  pf_r = mf_r.fem_of_element(cv);
1418  pf_m = mf_mult.fem_of_element(cv);
1419 
1420  if (!pf_m->is_lagrange() && !warning_msg1) {
1421  GMM_WARNING3("Dirichlet condition with non-lagrange multiplier fem. "
1422  "see the documentation about Dirichlet conditions.");
1423  warning_msg1 = true;
1424  }
1425 
1426  if (!(version & ASMDIR_SIMPLIFY)) continue;
1427 
1428  mesh_fem::ind_dof_face_ct pf_u_ct
1429  = mf_u.ind_basic_dof_of_face_of_element(cv, f);
1430  mesh_fem::ind_dof_face_ct pf_r_ct
1431  = mf_r.ind_basic_dof_of_face_of_element(cv, f);
1432  mesh_fem::ind_dof_face_ct pf_m_ct
1433  = mf_mult.ind_basic_dof_of_face_of_element(cv, f);
1434 
1435  size_type pf_u_nbdf = pf_u_ct.size();
1436  size_type pf_m_nbdf = pf_m_ct.size();
1437  size_type pf_u_nbdf_loc = pf_u->structure(cv)->nb_points_of_face(f);
1438  size_type pf_m_nbdf_loc = pf_m->structure(cv)->nb_points_of_face(f);
1439  // size_type pf_r_nbdf_loc = pf_r->structure(cv)->nb_points_of_face(f);
1440 
1441  if (pf_u_nbdf < pf_m_nbdf && !warning_msg2) {
1442  GMM_WARNING2("Dirichlet condition with a too rich multiplier fem. "
1443  "see the documentation about Dirichlet conditions.");
1444  warning_msg2 = true;
1445  }
1446 
1447  if (pf_u != pf_r || pf_u_nbdf != pf_m_nbdf ||
1448  ((pf_u != pf_r) && (pf_u_nbdf_loc != pf_m_nbdf_loc))) {
1449  for (size_type i = 0; i < pf_m_nbdf; ++i)
1450  nonsimplifiable_dofs.add(pf_m_ct[i]);
1451  continue;
1452  }
1453 
1454  for (size_type i = 0; i < pf_m_nbdf; ++i) {
1455  simplifiable_dofs.add(pf_m_ct[i]);
1456  simplifiable_indices[pf_m_ct[i]] = pf_u_ct[i];
1457  }
1458 
1459  if (!(version & ASMDIR_BUILDR)) continue;
1460 
1461  if (pf_u == pf_r) { // simplest simplification.
1462  size_type Qratio = mf_u.get_qdim() / mf_r.get_qdim();
1463  for (size_type i = 0; i < pf_m_nbdf; ++i) {
1464  simplifiable_values[pf_m_ct[i]]
1465  = r_data[pf_r_ct[i/Qratio]*Qratio+(i%Qratio)];
1466  }
1467  }
1468  }
1469 
1470  if (version & ASMDIR_SIMPLIFY) {
1471  if (simplifiable_dofs.card() > 0)
1472  { GMM_TRACE3("Simplification of the Dirichlet condition"); }
1473  else
1474  GMM_TRACE3("Sorry, no simplification of the Dirichlet condition");
1475  if (nonsimplifiable_dofs.card() > 0 && simplifiable_dofs.card() > 0)
1476  GMM_WARNING3("Partial simplification of the Dirichlet condition");
1477 
1478  for (dal::bv_visitor i(simplifiable_dofs); !i.finished(); ++i)
1479  if (!(nonsimplifiable_dofs[i])) {
1480  if (version & ASMDIR_BUILDH) { /* "erase" the row i */
1481  const mesh::ind_cv_ct &cv_ct = mf_mult.convex_to_basic_dof(i);
1482  for (size_type j = 0; j < cv_ct.size(); ++j) {
1483  size_type cv = cv_ct[j];
1484  for (size_type k=0; k < mf_u.nb_basic_dof_of_element(cv); ++k)
1485  H(i, mf_u.ind_basic_dof_of_element(cv)[k]) = value_type(0);
1486  }
1487  H(i, simplifiable_indices[i]) = value_type(1);
1488  }
1489  if (version & ASMDIR_BUILDR) R[i] = simplifiable_values[i];
1490  }
1491  }
1492  }
1493 
1494  template<typename MATRM, typename VECT1, typename VECT2>
1495  void assembling_Dirichlet_condition
1496  (MATRM &RM, VECT1 &B, const getfem::mesh_fem &mf, size_type boundary,
1497  const VECT2 &F) {
1498  // Works only for Lagrange dofs.
1499  size_type Q=mf.get_qdim();
1500  GMM_ASSERT1(!(mf.is_reduced()), "This function is not adapted to "
1501  "reduced finite element methods");
1502  dal::bit_vector nndof = mf.basic_dof_on_region(boundary);
1503  getfem::pfem pf1;
1504  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
1505  pf1 = mf.fem_of_element(cv);
1506  getfem::pdof_description ldof = getfem::lagrange_dof(pf1->dim());
1507  size_type nbd = pf1->nb_dof(cv);
1508  for (size_type i = 0; i < nbd; i++) {
1509  size_type dof1 = mf.ind_basic_dof_of_element(cv)[i*Q];
1510  if (nndof.is_in(dof1) && pf1->dof_types()[i] == ldof) {
1511  // cout << "dof : " << i << endl;
1512  for (size_type j = 0; j < nbd; j++) {
1513  size_type dof2 = mf.ind_basic_dof_of_element(cv)[j*Q];
1514  for (size_type k = 0; k < Q; ++k)
1515  for (size_type l = 0; l < Q; ++l) {
1516  if (!(nndof.is_in(dof2)) &&
1517  getfem::dof_compatibility(pf1->dof_types()[j],
1518  getfem::lagrange_dof(pf1->dim())))
1519  B[dof2+k] -= RM(dof2+k, dof1+l) * F[dof1+l];
1520  RM(dof2+k, dof1+l) = RM(dof1+l, dof2+k) = 0;
1521  }
1522  }
1523  for (size_type k = 0; k < Q; ++k)
1524  { RM(dof1+k, dof1+k) = 1; B[dof1+k] = F[dof1+k]; }
1525  }
1526  }
1527  }
1528  }
1529 
1530  /**
1531  Assembly of generalized Dirichlet constraints h(x)u(x) = r(x),
1532  where h is a QxQ matrix field (Q == mf_u.get_qdim()), outputs a
1533  (under-determined) linear system HU=R.
1534 
1535  size(h_data) = Q^2 * nb_dof(mf_rh);
1536  size(r_data) = Q * nb_dof(mf_rh);
1537 
1538  This function tries hard to make H diagonal or mostly diagonal:
1539  this function is able to "simplify" the dirichlet constraints (see below)
1540  version = |ASMDIR_BUILDH : build H
1541  |ASMDIR_BUILDR : build R
1542  |ASMDIR_SIMPLIFY : simplify
1543  |ASMDIR_BUILDALL : do everything.
1544 
1545  @ingroup asm
1546  */
1547 
1548  template<typename MAT, typename VECT1, typename VECT2, typename VECT3>
1549  void asm_generalized_dirichlet_constraints
1550  (MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u,
1551  const mesh_fem &mf_h, const mesh_fem &mf_r, const VECT2 &h_data,
1552  const VECT3 &r_data, const mesh_region &region,
1553  int version = ASMDIR_BUILDALL) {
1554  pfem pf_u, pf_rh;
1555 
1556  if ((version & ASMDIR_SIMPLIFY) &&
1557  (mf_u.is_reduced() || mf_h.is_reduced() || mf_r.is_reduced())) {
1558  GMM_WARNING1("Sorry, no simplification for reduced fems");
1559  version = (version & ASMDIR_BUILDR);
1560  }
1561 
1562  region.from_mesh(mim.linked_mesh()).error_if_not_faces();
1563  GMM_ASSERT1(mf_h.get_qdim() == 1 && mf_r.get_qdim() == 1,
1564  "invalid data mesh fem (Qdim=1 required)");
1565  if (version & ASMDIR_BUILDH) {
1566  asm_qu_term(H, mim, mf_u, mf_h, h_data, region);
1567  std::vector<size_type> ind(0);
1568  dal::bit_vector bdof = mf_u.basic_dof_on_region(region);
1569  // gmm::clean(H, 1E-15 * gmm::mat_maxnorm(H));
1570  for (size_type i = 0; i < mf_u.nb_dof(); ++i)
1571  if (!(bdof[i])) ind.push_back(i);
1572  gmm::clear(gmm::sub_matrix(H, gmm::sub_index(ind)));
1573  }
1574  if (version & ASMDIR_BUILDR)
1575  asm_source_term(R, mim, mf_u, mf_r, r_data, region);
1576  if (!(version & ASMDIR_SIMPLIFY)) return;
1577 
1578  /* step 2 : simplification of simple dirichlet conditions */
1579  if (&mf_r == &mf_h) {
1580  for (mr_visitor v(region); !v.finished(); v.next()) {
1581  size_type cv = v.cv();
1582  short_type f = v.f();
1583 
1584  GMM_ASSERT1(mf_u.convex_index().is_in(cv) &&
1585  mf_r.convex_index().is_in(cv),
1586  "attempt to impose a dirichlet "
1587  "condition on a convex with no FEM!");
1588 
1589  if (f >= mf_u.linked_mesh().structure_of_convex(cv)->nb_faces())
1590  continue;
1591  pf_u = mf_u.fem_of_element(cv);
1592  pf_rh = mf_r.fem_of_element(cv);
1593  /* don't try anything with vector elements */
1594  if (mf_u.fem_of_element(cv)->target_dim() != 1) continue;
1595  bgeot::pconvex_structure cvs_u = pf_u->structure(cv);
1596  bgeot::pconvex_structure cvs_rh = pf_rh->structure(cv);
1597  for (size_type i = 0; i < cvs_u->nb_points_of_face(f); ++i) {
1598 
1599  size_type Q = mf_u.get_qdim(); // pf_u->target_dim() (==1)
1600 
1601  size_type ind_u = cvs_u->ind_points_of_face(f)[i];
1602  pdof_description tdof_u = pf_u->dof_types()[ind_u];
1603 
1604  for (size_type j = 0; j < cvs_rh->nb_points_of_face(f); ++j) {
1605  size_type ind_rh = cvs_rh->ind_points_of_face(f)[j];
1606  pdof_description tdof_rh = pf_rh->dof_types()[ind_rh];
1607  /*
1608  same kind of dof and same location of dof ?
1609  => then the previous was not useful for this dofs (introducing
1610  a mass matrix which is not diagonal in the constraints matrix)
1611  -> the constraint is simplified:
1612  we replace \int{(H_j.psi_j)*phi_i}=\int{R_j.psi_j} (sum over j)
1613  with H_j*phi_i = R_j
1614  --> Le principe peut être faux : non identique à la projection
1615  L^2 et peut entrer en conccurence avec les autres ddl -> a revoir
1616  */
1617  if (tdof_u == tdof_rh &&
1618  gmm::vect_dist2_sqr((*(pf_u->node_tab(cv)))[ind_u],
1619  (*(pf_rh->node_tab(cv)))[ind_rh])
1620  < 1.0E-14) {
1621  /* the dof might be "duplicated" */
1622  for (size_type q = 0; q < Q; ++q) {
1623  size_type dof_u = mf_u.ind_basic_dof_of_element(cv)[ind_u*Q+q];
1624 
1625  /* "erase" the row */
1626  if (version & ASMDIR_BUILDH)
1627  for (size_type k=0; k < mf_u.nb_basic_dof_of_element(cv); ++k)
1628  H(dof_u, mf_u.ind_basic_dof_of_element(cv)[k]) = 0.0;
1629 
1630  size_type dof_rh = mf_r.ind_basic_dof_of_element(cv)[ind_rh];
1631  /* set the "simplified" row */
1632  if (version & ASMDIR_BUILDH)
1633  for (unsigned jj=0; jj < Q; jj++) {
1634  size_type dof_u2
1635  = mf_u.ind_basic_dof_of_element(cv)[ind_u*Q+jj];
1636  H(dof_u, dof_u2) = h_data[(jj*Q+q) + Q*Q*(dof_rh)];
1637  }
1638  if (version & ASMDIR_BUILDR) R[dof_u] = r_data[dof_rh*Q+q];
1639  }
1640  }
1641  }
1642  }
1643  }
1644  }
1645  }
1646 
1647  /**
1648  Build an orthogonal basis of the kernel of H in NS, gives the
1649  solution of minimal norm of H*U = R in U0 and return the
1650  dimension of the kernel. The function is based on a
1651  Gramm-Schmidt algorithm.
1652 
1653  @ingroup asm
1654  */
1655  template<typename MAT1, typename MAT2, typename VECT1, typename VECT2>
1656  size_type Dirichlet_nullspace(const MAT1 &H, MAT2 &NS,
1657  const VECT1 &R, VECT2 &U0) {
1658 
1659  typedef typename gmm::linalg_traits<MAT1>::value_type T;
1660  typedef typename gmm::number_traits<T>::magnitude_type MAGT;
1661  typedef typename gmm::temporary_vector<MAT1>::vector_type TEMP_VECT;
1662  MAGT tol = gmm::default_tol(MAGT());
1663  MAGT norminfH = gmm::mat_maxnorm(H);
1664  size_type nbd = gmm::mat_ncols(H), nbase = 0, nbr = gmm::mat_nrows(H);
1665  TEMP_VECT aux(nbr), e(nbd), f(nbd);
1666  dal::dynamic_array<TEMP_VECT> base_img;
1667  dal::dynamic_array<TEMP_VECT> base_img_inv;
1668  size_type nb_bimg = 0;
1669  gmm::clear(NS);
1670 
1671  if (!(gmm::is_col_matrix(H)))
1672  GMM_WARNING2("Dirichlet_nullspace inefficient for a row matrix H");
1673  // First, detection of null columns of H, and already orthogonals
1674  // vectors of the image of H.
1675  dal::bit_vector nn;
1676  for (size_type i = 0; i < nbd; ++i) {
1677  gmm::clear(e); e[i] = T(1);
1678  gmm::mult(H, e, aux);
1679  MAGT n = gmm::vect_norm2(aux);
1680 
1681  if (n < norminfH*tol*1000) {
1682  NS(i, nbase++) = T(1); nn[i] = true;
1683  }
1684  else {
1685  bool good = true;
1686  for (size_type j = 0; j < nb_bimg; ++j)
1687  if (gmm::abs(gmm::vect_sp(aux, base_img[j])) > MAGT(0))
1688  { good = false; break; }
1689  if (good) {
1690  gmm::copy(e, f);
1691  gmm::scale(f, T(MAGT(1)/n)); gmm::scale(aux, T(MAGT(1)/n));
1692  base_img_inv[nb_bimg] = TEMP_VECT(nbd);
1693  gmm::copy(f, base_img_inv[nb_bimg]);
1694  gmm::clean(aux, gmm::vect_norminf(aux)*tol);
1695  base_img[nb_bimg] = TEMP_VECT(nbr);
1696  gmm::copy(aux, base_img[nb_bimg++]);
1697  nn[i] = true;
1698  }
1699  }
1700  }
1701  size_type nb_triv_base = nbase;
1702 
1703  for (size_type i = 0; i < nbd; ++i) {
1704  if (!(nn[i])) {
1705  gmm::clear(e); e[i] = T(1); gmm::clear(f); f[i] = T(1);
1706  gmm::mult(H, e, aux);
1707  for (size_type j = 0; j < nb_bimg; ++j) {
1708  T c = gmm::vect_sp(aux, base_img[j]);
1709  // if (gmm::abs(c > 1.0E-6) { // à scaler sur l'ensemble de H ...
1710  if (c != T(0)) {
1711  gmm::add(gmm::scaled(base_img[j], -c), aux);
1712  gmm::add(gmm::scaled(base_img_inv[j], -c), f);
1713  }
1714  }
1715  if (gmm::vect_norm2(aux) < norminfH*tol*MAGT(10000)) {
1716  gmm::copy(f, gmm::mat_col(NS, nbase++));
1717  }
1718  else {
1719  MAGT n = gmm::vect_norm2(aux);
1720  gmm::scale(f, T(MAGT(1)/n)); gmm::scale(aux, T(MAGT(1)/n));
1721  gmm::clean(f, tol*gmm::vect_norminf(f));
1722  gmm::clean(aux, tol*gmm::vect_norminf(aux));
1723  base_img_inv[nb_bimg] = TEMP_VECT(nbd);
1724  gmm::copy(f, base_img_inv[nb_bimg]);
1725  base_img[nb_bimg] = TEMP_VECT(nbr);
1726  gmm::copy(aux, base_img[nb_bimg++]);
1727  }
1728  }
1729  }
1730  // Compute a solution in U0
1731  gmm::clear(U0);
1732  for (size_type i = 0; i < nb_bimg; ++i) {
1733  T c = gmm::vect_sp(base_img[i], R);
1734  gmm::add(gmm::scaled(base_img_inv[i], c), U0);
1735  }
1736  // Orthogonalisation of the basis of the kernel of H.
1737  for (size_type i = nb_triv_base; i < nbase; ++i) {
1738  for (size_type j = nb_triv_base; j < i; ++j) {
1739  T c = gmm::vect_sp(gmm::mat_col(NS,i), gmm::mat_col(NS,j));
1740  if (c != T(0))
1741  gmm::add(gmm::scaled(gmm::mat_col(NS,j), -c), gmm::mat_col(NS,i));
1742  }
1743  gmm::scale(gmm::mat_col(NS,i),
1744  T(1) / gmm::vect_norm2(gmm::mat_col(NS,i)));
1745  }
1746  // projection of U0 on the orthogonal to the kernel.
1747  for (size_type j = nb_triv_base; j < nbase; ++j) {
1748  T c = gmm::vect_sp(gmm::mat_col(NS,j), U0);
1749  if (c != T(0))
1750  gmm::add(gmm::scaled(gmm::mat_col(NS,j), -c), U0);
1751  }
1752  // Test ...
1753  gmm::mult(H, U0, gmm::scaled(R, T(-1)), aux);
1754  if (gmm::vect_norm2(aux) > gmm::vect_norm2(U0)*tol*MAGT(10000))
1755  GMM_WARNING2("Dirichlet condition not well inverted: residu="
1756  << gmm::vect_norm2(aux));
1757 
1758  return nbase;
1759  }
1760 
1761 } /* end of namespace getfem. */
1762 
1763 
1764 #endif /* GETFEM_ASSEMBLING_H__ */
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
dal::dynamic_array
Dynamic Array.
Definition: dal_basic.h:196
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:152
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:454
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:316
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:566
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_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:560
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::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:472
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!...
Definition: getfem_mesh_fem.h:448
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:195
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:498
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:199
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::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
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::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:237
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:55
getfem::mesh_region::from_mesh
const mesh_region & from_mesh(const mesh &m) const
For regions which have been built with just a number 'id', from_mesh(m) sets the current region to 'm...
Definition: getfem_mesh_region.cc:72
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:142
getfem_assembling_tensors.h
Generic assembly implementation.
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1277
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_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_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:1016
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_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::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:1110
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:1152
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:1195
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:1656
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:974
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_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_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::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_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:1344
getfem::asm_lumped_mass_matrix_for_first_order_param
void asm_lumped_mass_matrix_for_first_order_param(MAT &M, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_data, const VECT &F, const mesh_region &rg=mesh_region::all_convexes())
lumped mass matrix assembly with an additional parameter (on the whole mesh or on the specified bound...
Definition: getfem_assembling.h:866
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::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:1274
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
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
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:1140
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:1063
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.
Definition: getfem_assembling.h:1373
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::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::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_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:1074
getfem::asm_lumped_mass_matrix_for_first_order
void asm_lumped_mass_matrix_for_first_order(const MAT &M, const mesh_im &mim, const mesh_fem &mf1, const mesh_region &rg=mesh_region::all_convexes())
lumped mass matrix assembly (on the whole mesh or on the specified boundary)
Definition: getfem_assembling.h:853
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_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:1550
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
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:944
getfem::dof_compatibility
bool dof_compatibility(pdof_description, pdof_description)
Says if the two dofs can be identified.
Definition: getfem_fem.cc:619
getfem::lagrange_dof
pdof_description lagrange_dof(dim_type d)
Description of a unique dof of lagrange type (value at the node).
Definition: getfem_fem.cc:391
getfem::pdof_description
dof_description * pdof_description
Type representing a pointer on a dof_description.
Definition: getfem_fem.h:160
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:73
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
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::asm_lumped_mass_matrix_for_first_order_from_consistent
void asm_lumped_mass_matrix_for_first_order_from_consistent(const MAT &M)
lumped mass matrix assembly from consistent mass matrix
Definition: getfem_assembling.h:833
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:1229
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:1217
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:1243
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:1165
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:1206
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:1257