1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributor(s):
7 // Eric Bechet
8 //
9
10 #ifndef TERMS_H
11 #define TERMS_H
12
13 #include "SVector3.h"
14 #include "STensor3.h"
15 #include "STensor43.h"
16 #include "Numeric.h"
17 #include "functionSpace.h"
18 #include "groupOfElements.h"
19 #include "materialLaw.h"
20 #include <vector>
21 #include <iterator>
22
23 template <class T2> class ScalarTermBase;
24
25 template <class T2> class LinearTermBase;
26 template <class T2> class PlusTerm;
27
28 class BilinearTermBase;
29
dot(const double & a,const double & b)30 inline double dot(const double &a, const double &b) { return a * b; }
31
delta(int i,int j)32 inline int delta(int i, int j)
33 {
34 if(i == j)
35 return 1;
36 else
37 return 0;
38 }
39
setDirection(double & a,const double & b)40 inline void setDirection(double &a, const double &b) { a = b; }
setDirection(SVector3 & a,const SVector3 & b)41 inline void setDirection(SVector3 &a, const SVector3 &b)
42 {
43 for(int i = 0; i < 3; i++) a(i) = b(i);
44 }
45
46 template <class T2 = double> class ScalarTermBase {
47 public:
~ScalarTermBase()48 virtual ~ScalarTermBase() {}
49 virtual void get(MElement *ele, int npts, IntPt *GP, T2 &val) const = 0;
50 virtual void get(MElement *ele, int npts, IntPt *GP,
51 std::vector<T2> &vval) const = 0;
52 virtual ScalarTermBase<T2> *clone() const = 0;
53 };
54
55 template <class T2 = double> class ScalarTerm : public ScalarTermBase<T2> {
56 public:
~ScalarTerm()57 virtual ~ScalarTerm() {}
58 };
59
60 template <class T2 = double> class ScalarTermConstant : public ScalarTerm<T2> {
61 protected:
62 T2 cst;
63
64 public:
cst(val_)65 ScalarTermConstant(T2 val_ = T2()) : cst(val_) {}
~ScalarTermConstant()66 virtual ~ScalarTermConstant() {}
67 virtual void get(MElement *ele, int npts, IntPt *GP, T2 &val) const;
68 virtual void get(MElement *ele, int npts, IntPt *GP,
69 std::vector<T2> &vval) const;
70 virtual void get(MVertex *ver, T2 &val) const;
clone()71 virtual ScalarTermBase<T2> *clone() const
72 {
73 return new ScalarTermConstant<T2>(cst);
74 }
75 };
76
77 class BilinearTermToScalarTerm : public ScalarTerm<double> {
78 protected:
79 BilinearTermBase &bilterm;
80
81 public:
BilinearTermToScalarTerm(BilinearTermBase & bilterm_)82 BilinearTermToScalarTerm(BilinearTermBase &bilterm_) : bilterm(bilterm_) {}
~BilinearTermToScalarTerm()83 virtual ~BilinearTermToScalarTerm() {}
84 virtual void get(MElement *ele, int npts, IntPt *GP, double &val) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<double> & vval)85 virtual void get(MElement *ele, int npts, IntPt *GP,
86 std::vector<double> &vval) const {};
clone()87 virtual ScalarTermBase<double> *clone() const
88 {
89 return new BilinearTermToScalarTerm(bilterm);
90 }
91 };
92
93 class BilinearTermBase {
94 public:
~BilinearTermBase()95 virtual ~BilinearTermBase() {}
96 virtual void get(MElement *ele, int npts, IntPt *GP,
97 fullMatrix<double> &m) const;
98 virtual void get(MElement *ele, int npts, IntPt *GP,
99 std::vector<fullMatrix<double> > &mv) const = 0;
100 virtual BilinearTermBase *clone() const = 0;
101 };
102
103 template <class T2> class BilinearTermContract : public BilinearTermBase {
104 protected:
105 const LinearTermBase<T2> *a;
106 const LinearTermBase<T2> *b;
107
108 public:
BilinearTermContract(const LinearTermBase<T2> & a_,const LinearTermBase<T2> & b_)109 BilinearTermContract(const LinearTermBase<T2> &a_,
110 const LinearTermBase<T2> &b_)
111 : a(a_.clone()), b(b_.clone())
112 {
113 }
~BilinearTermContract()114 virtual ~BilinearTermContract()
115 {
116 delete a;
117 delete b;
118 }
119 virtual void get(MElement *ele, int npts, IntPt *GP,
120 fullMatrix<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullMatrix<double>> & mv)121 virtual void get(MElement *ele, int npts, IntPt *GP,
122 std::vector<fullMatrix<double> > &mv) const {};
clone()123 virtual BilinearTermBase *clone() const
124 {
125 return new BilinearTermContract<T2>(*a, *b);
126 }
127 };
128
129 template <class T2>
130 class BilinearTermContractWithLaw : public BilinearTermContract<T2> {
131 public:
132 protected:
133 const ScalarTermBase<typename TensorialTraits<T2>::TensProdType> *c;
134
135 public:
BilinearTermContractWithLaw(const LinearTermBase<T2> & a_,const ScalarTermBase<typename TensorialTraits<T2>::TensProdType> & c_,const LinearTermBase<T2> & b_)136 BilinearTermContractWithLaw(
137 const LinearTermBase<T2> &a_,
138 const ScalarTermBase<typename TensorialTraits<T2>::TensProdType> &c_,
139 const LinearTermBase<T2> &b_)
140 : BilinearTermContract<T2>(a_, b_), c(c_.clone())
141 {
142 }
~BilinearTermContractWithLaw()143 virtual ~BilinearTermContractWithLaw() { delete c; }
144 virtual void get(MElement *ele, int npts, IntPt *GP,
145 fullMatrix<double> &m) const;
146 virtual void get(MElement *ele, int npts, IntPt *GP,
147 std::vector<fullMatrix<double> > &mv) const;
clone()148 virtual BilinearTermBase *clone() const
149 {
150 return new BilinearTermContractWithLaw<T2>(*BilinearTermContract<T2>::a, *c,
151 *BilinearTermContract<T2>::b);
152 }
153 };
154
155 template <class T2>
156 BilinearTermContract<T2> operator|(const LinearTermBase<T2> &L1,
157 const LinearTermBase<T2> &L2);
158
159 template <class T1, class T2> class BilinearTerm : public BilinearTermBase {
160 protected:
161 FunctionSpace<T1> &space1;
162 FunctionSpace<T2> &space2;
163
164 public:
BilinearTerm(FunctionSpace<T1> & space1_,FunctionSpace<T2> & space2_)165 BilinearTerm(FunctionSpace<T1> &space1_, FunctionSpace<T2> &space2_)
166 : space1(space1_), space2(space2_)
167 {
168 }
~BilinearTerm()169 virtual ~BilinearTerm() {}
170 };
171
172 template <class T2 = double> class LinearTermBase {
173 public:
~LinearTermBase()174 virtual ~LinearTermBase() {}
175 virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const;
176 virtual void get(MElement *ele, int npts, IntPt *GP,
177 std::vector<fullVector<T2> > &vv) const = 0;
178 virtual LinearTermBase<T2> *clone() const = 0;
179 PlusTerm<T2> operator+(const LinearTermBase<T2> &other);
180 };
181
182 template <class T2> class PlusTerm : public LinearTermBase<T2> {
183 const LinearTermBase<T2> *a;
184 const LinearTermBase<T2> *b;
185
186 public:
PlusTerm(const LinearTermBase<T2> & a_,const LinearTermBase<T2> & b_)187 PlusTerm(const LinearTermBase<T2> &a_, const LinearTermBase<T2> &b_)
188 : a(a_.clone()), b(b_.clone())
189 {
190 }
~PlusTerm()191 virtual ~PlusTerm()
192 {
193 delete a;
194 delete b;
195 }
196 virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullVector<T2>> & vv)197 virtual void get(MElement *ele, int npts, IntPt *GP,
198 std::vector<fullVector<T2> > &vv) const {};
clone()199 virtual LinearTermBase<T2> *clone() const { return new PlusTerm<T2>(*a, *b); }
200 };
201
202 template <class T1, class T2 = double>
203 class LinearTerm : public LinearTermBase<T2> {
204 protected:
205 FunctionSpace<T1> &space1;
206
207 public:
LinearTerm(FunctionSpace<T1> & space1_)208 LinearTerm(FunctionSpace<T1> &space1_) : space1(space1_) {}
~LinearTerm()209 virtual ~LinearTerm() {}
210 };
211
212 template <class T1>
213 class GradTerm : public LinearTerm<T1, typename TensorialTraits<T1>::GradType> {
214 public:
GradTerm(FunctionSpace<T1> & space1_)215 GradTerm(FunctionSpace<T1> &space1_)
216 : LinearTerm<T1, typename TensorialTraits<T1>::GradType>(space1_)
217 {
218 }
219 virtual void
get(MElement * ele,int npts,IntPt * GP,fullVector<typename TensorialTraits<T1>::GradType> & vec)220 get(MElement *ele, int npts, IntPt *GP,
221 fullVector<typename TensorialTraits<T1>::GradType> &vec) const
222 {
223 LinearTermBase<typename TensorialTraits<T1>::GradType>::get(ele, npts, GP,
224 vec);
225 }
226 virtual void
227 get(MElement *ele, int npts, IntPt *GP,
228 std::vector<fullVector<typename TensorialTraits<T1>::GradType> > &vvec)
229 const;
clone()230 virtual LinearTermBase<typename TensorialTraits<T1>::GradType> *clone() const
231 {
232 return new GradTerm<T1>(
233 LinearTerm<T1, typename TensorialTraits<T1>::GradType>::space1);
234 }
~GradTerm()235 virtual ~GradTerm() {}
236 };
237
238 template <class T1, class T2> class LaplaceTerm : public BilinearTerm<T1, T2> {
239 public:
LaplaceTerm(FunctionSpace<T1> & space1_,FunctionSpace<T2> & space2_)240 LaplaceTerm(FunctionSpace<T1> &space1_, FunctionSpace<T2> &space2_)
241 : BilinearTerm<T1, T2>(space1_, space2_)
242 {
243 }
~LaplaceTerm()244 virtual ~LaplaceTerm() {}
get(MElement * ele,int npts,IntPt * GP,fullMatrix<double> & m)245 virtual void get(MElement *ele, int npts, IntPt *GP,
246 fullMatrix<double> &m) const
247 {
248 Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
249 }
get(MElement * ele,int npts,IntPt * GP,std::vector<fullMatrix<double>> & vm)250 virtual void get(MElement *ele, int npts, IntPt *GP,
251 std::vector<fullMatrix<double> > &vm) const
252 {
253 Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
254 }
get(MVertex * ver,fullMatrix<double> & m)255 virtual void get(MVertex *ver, fullMatrix<double> &m)
256 {
257 Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
258 }
clone()259 virtual BilinearTermBase *clone() const
260 {
261 return new LaplaceTerm<T1, T2>(BilinearTerm<T1, T2>::space1,
262 BilinearTerm<T1, T2>::space2);
263 }
264 }; // class
265
266 template <class T1>
267 class LaplaceTerm<T1, T1>
268 : public BilinearTerm<T1, T1> // symmetric
269 {
270 protected:
271 double diffusivity;
272
273 public:
274 LaplaceTerm(FunctionSpace<T1> &space1_, double diff = 1)
275 : BilinearTerm<T1, T1>(space1_, space1_), diffusivity(diff)
276 {
277 }
~LaplaceTerm()278 virtual ~LaplaceTerm() {}
279 virtual void get(MElement *ele, int npts, IntPt *GP,
280 fullMatrix<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullMatrix<double>> & vm)281 virtual void get(MElement *ele, int npts, IntPt *GP,
282 std::vector<fullMatrix<double> > &vm) const {};
clone()283 virtual BilinearTermBase *clone() const
284 {
285 return new LaplaceTerm<T1, T1>(BilinearTerm<T1, T1>::space1, diffusivity);
286 }
287 }; // class
288
289 class IsotropicElasticTerm : public BilinearTerm<SVector3, SVector3> {
290 protected:
291 double E, nu;
292 bool sym;
293 fullMatrix<double> H; /* =
294 { {C11, C12, C12, 0, 0, 0},
295 {C12, C11, C12, 0, 0, 0},
296 {C12, C12, C11, 0, 0, 0},
297 { 0, 0, 0, C44, 0, 0},
298 { 0, 0, 0, 0, C44, 0},
299 { 0, 0, 0, 0, 0, C44} };*/
300 public:
301 IsotropicElasticTerm(FunctionSpace<SVector3> &space1_,
302 FunctionSpace<SVector3> &space2_, double E_, double nu_);
303 IsotropicElasticTerm(FunctionSpace<SVector3> &space1_, double E_, double nu_);
~IsotropicElasticTerm()304 virtual ~IsotropicElasticTerm() {}
305 virtual void get(MElement *ele, int npts, IntPt *GP,
306 fullMatrix<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullMatrix<double>> & vm)307 virtual void get(MElement *ele, int npts, IntPt *GP,
308 std::vector<fullMatrix<double> > &vm) const {};
clone()309 virtual BilinearTermBase *clone() const
310 {
311 return new IsotropicElasticTerm(BilinearTerm<SVector3, SVector3>::space1,
312 BilinearTerm<SVector3, SVector3>::space2, E,
313 nu);
314 }
315 }; // class
316
317 template <class T1> class LoadTerm : public LinearTerm<T1> {
318 protected:
319 simpleFunction<typename TensorialTraits<T1>::ValType> *Load;
320
321 public:
LoadTerm(FunctionSpace<T1> & space1_,simpleFunction<typename TensorialTraits<T1>::ValType> * Load_)322 LoadTerm(FunctionSpace<T1> &space1_,
323 simpleFunction<typename TensorialTraits<T1>::ValType> *Load_)
324 : LinearTerm<T1>(space1_), Load(Load_)
325 {
326 }
~LoadTerm()327 virtual ~LoadTerm() {}
clone()328 virtual LinearTermBase<double> *clone() const
329 {
330 return new LoadTerm<T1>(LinearTerm<T1>::space1, Load);
331 }
332 virtual void get(MElement *ele, int npts, IntPt *GP,
333 fullVector<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullVector<double>> & vv)334 virtual void get(MElement *ele, int npts, IntPt *GP,
335 std::vector<fullVector<double> > &vv) const {};
336 };
337
338 template <class T1>
339 class LagrangeMultiplierTerm : public BilinearTerm<T1, double> {
340 T1 _d;
341
342 public:
LagrangeMultiplierTerm(FunctionSpace<T1> & space1_,FunctionSpace<double> & space2_,const T1 & d)343 LagrangeMultiplierTerm(FunctionSpace<T1> &space1_,
344 FunctionSpace<double> &space2_, const T1 &d)
345 : BilinearTerm<T1, double>(space1_, space2_)
346 {
347 setDirection(_d, d);
348 }
~LagrangeMultiplierTerm()349 virtual ~LagrangeMultiplierTerm() {}
350 virtual void get(MElement *ele, int npts, IntPt *GP,
351 fullMatrix<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullMatrix<double>> & vm)352 virtual void get(MElement *ele, int npts, IntPt *GP,
353 std::vector<fullMatrix<double> > &vm) const {};
clone()354 virtual BilinearTermBase *clone() const
355 {
356 return new LagrangeMultiplierTerm(BilinearTerm<T1, double>::space1,
357 BilinearTerm<T1, double>::space2, _d);
358 }
359 };
360
361 class LagMultTerm : public BilinearTerm<SVector3, SVector3> {
362 private:
363 double _eqfac;
364
365 public:
366 LagMultTerm(FunctionSpace<SVector3> &space1_,
367 FunctionSpace<SVector3> &space2_, double eqfac = 1.0)
368 : BilinearTerm<SVector3, SVector3>(space1_, space2_), _eqfac(eqfac)
369 {
370 }
~LagMultTerm()371 virtual ~LagMultTerm() {}
372 virtual void get(MElement *ele, int npts, IntPt *GP,
373 fullMatrix<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullMatrix<double>> & vm)374 virtual void get(MElement *ele, int npts, IntPt *GP,
375 std::vector<fullMatrix<double> > &vm) const {};
clone()376 virtual BilinearTermBase *clone() const
377 {
378 return new LagMultTerm(BilinearTerm<SVector3, SVector3>::space1,
379 BilinearTerm<SVector3, SVector3>::space2, _eqfac);
380 }
381 };
382
383 template <class T1> class LoadTermOnBorder : public LinearTerm<T1> {
384 private:
385 double _eqfac;
386 simpleFunction<typename TensorialTraits<T1>::ValType> *Load;
387
388 public:
389 LoadTermOnBorder(FunctionSpace<T1> &space1_,
390 simpleFunction<typename TensorialTraits<T1>::ValType> *Load_,
391 double eqfac = 1.0)
392 : LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_)
393 {
394 }
~LoadTermOnBorder()395 virtual ~LoadTermOnBorder() {}
clone()396 virtual LinearTermBase<double> *clone() const
397 {
398 return new LoadTermOnBorder<T1>(LinearTerm<T1>::space1, Load, _eqfac);
399 }
400 virtual void get(MElement *ele, int npts, IntPt *GP,
401 fullVector<double> &m) const;
get(MElement * ele,int npts,IntPt * GP,std::vector<fullVector<double>> & vv)402 virtual void get(MElement *ele, int npts, IntPt *GP,
403 std::vector<fullVector<double> > &vv) const {};
404 };
405
406 #include "terms.hpp"
407
408 #endif
409