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