1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file OperatorOnUnknown.hpp
19 \author E. Lunéville
20 \since 02 mar 2012
21 \date 9 may 2012
22
23 \brief Definition of the xlifepp::OperatorOnUnknown class
24
25 class xlifepp::OperatorOnUnknown stores the action of a differential operator on an unknown or test function
26 the general following form is accepted :
27 \verbatim
28 left_operand left_op diff_op(unknown) right_op right_operand
29 \endverbatim
30
31 where
32 - diff_op is a differential operator of type defined in the DifferentialOperator class
33 - left_op and right_op are algebraic operators of type product (*), inner product (|), cross product (^) or contracted product (%)
34 - left_operand and right_operand are either Function object or Value object (see related classes)
35
36 All the attributes may be omitted, except unknown!
37
38 For instance, the following syntaxes are available
39 \verbatim
40 u (trivial form)
41 grad(u), div(u), ... (simple form)
42 A*grad(u) (A a constant scalar or matrix or a function returning a scalar or a Matrix)
43 F|grad(u) (F a constant vector or a function returning a vector)
44 (A*grad(u))|F (a combination of the two previous forms)
45 \endverbatim
46
47 The left and right data may be conjugate or transpose : (tran(A)*grad(u))|conj(F)
48 but not the unknown and not the whole expression or a part of it
49
50 It is possible to use predefined vectors such as the normal vector _n
51 As these forms are processed using overloaded operators, be cautious with the C++ operator priority.
52 In doubt, use paranthesis.
53
54 Note that some compatibility tests are performed but not all possible.
55 Generally, structure (scalar, vector or matrix) consistancy is checked but not the dimension consistancy!
56 Therefore, it is the responsability of the user to insure the consistancy of his expression
57
58 The class LcOperatorOnUnknown handles linear combination of OperatorOnUnknown's
59 and inherits from std::vector<OpuValPair> where OpuValPair is an alias of std::pair<OperatorOnUnknown*, complex_t>
60 */
61
62 #ifndef OPERATOR_ON_UNKNOWN_HPP
63 #define OPERATOR_ON_UNKNOWN_HPP
64
65 #include "config.h"
66 #include "DifferentialOperator.hpp"
67 #include "Operand.hpp"
68 #include "space.h"
69
70 namespace xlifepp
71 {
72
73 //forward declaration
74 class EssentialCondition;
75 class TermVector;
76
77 /*!
78 \class OperatorOnUnknown
79 describes the operations to process on unknown.
80 Unknown may be scalar or vector
81 */
82 class OperatorOnUnknown
83 {
84 protected :
85 const Unknown* u_p; //!< unknown involved in operator
86 bool conjugateUnknown_; //!< true if the unknown has to be conjugated
87 DifferentialOperator* difOp_p; //!< differential operator involved in operator
88 Operand* leftOperand_p; //!< object before the diff operator
89 Operand* rightOperand_p; //!< object after the diff operator
90 std::vector<complex_t> coefs_; //!< coefficients for generalized operator (gradG, divG, curlG, epsilonG)
91 ValueType type_; //!< type of returned value (one among _real, _complex)
92 StrucType struct_; //!< structure of returned value (one among _scalar, _vector, _matrix)
93 bool leftPriority_; //!< priority order flag (true if left operand is applied first)
94 dimPair dimsRes_; //!< dimension of result, (0,0) if unset
95
96 void setReturnedType(const Unknown*, DiffOpType); //!< set returned type
97
98 public :
99 // constructors/destructor
100 OperatorOnUnknown(const Unknown* un = 0, DiffOpType = _id); //!< default constructor
101 OperatorOnUnknown(const Unknown&, DiffOpType = _id); //!< constructor with differential operator
102 OperatorOnUnknown(const Unknown&, DiffOpType, const std::vector<complex_t>&);//!< constructor with generalized differential operator
103 OperatorOnUnknown(const Unknown&, const Function&, AlgebraicOperator, bool); //!< useful constructor : (fun op u) or (u op fun)
104 OperatorOnUnknown(const Unknown&, const OperatorOnFunction&, AlgebraicOperator, bool); //!< useful constructor : (opfun op u) or (u op opfun)
105 OperatorOnUnknown(const Unknown&, const Value&, AlgebraicOperator, bool); //!< useful constructor : (val op u) or (u op val)
106 OperatorOnUnknown(const OperatorOnUnknown&); //!< copy constructor
107 ~OperatorOnUnknown(); //!< destructor
108 OperatorOnUnknown& operator=(const OperatorOnUnknown&); //!< assignment operator
109
110 void updateReturnedType(AlgebraicOperator, ValueType, StrucType, dimPair, bool); //!< update value and struct returned
111 void setStructure(); //!< set returned type, structure, dims of returned value
112
113 //accessors
unknown() const114 const Unknown* unknown() const {return u_p;} //!< return the unknown pointer
conjugateUnknown() const115 bool conjugateUnknown() const {return conjugateUnknown_;} //!< returns true if the unknown has to be conjugated
difOpType() const116 DiffOpType difOpType() const {return difOp_p->type();} //!< return type of differential operator
difOp_()117 DifferentialOperator*& difOp_() {return difOp_p;} //!< return pointer to differential operator
difOp() const118 DifferentialOperator& difOp() const {return *difOp_p;} //!< return reference to differential operator
valueType()119 ValueType& valueType() {return type_;} //!< return the operator returned type (non const)
valueType() const120 ValueType valueType()const {return type_;} //!< return the operator returned type
strucType() const121 StrucType strucType() const {return struct_;} //!< return the operator returned structure
dimsRes() const122 dimPair dimsRes() const {return dimsRes_;} //!< return the dimension of returned values
leftOperand()123 Operand*& leftOperand() {return leftOperand_p;} //!< return pointer to left operand
rightOperand()124 Operand*& rightOperand() {return rightOperand_p;} //!< return pointer to right operand
leftOperand() const125 const Operand* leftOperand() const //! return pointer to left operand
126 {return leftOperand_p;}
rightOperand() const127 const Operand* rightOperand() const //! return pointer to left operand
128 {return rightOperand_p;}
coefs()129 std::vector<complex_t>& coefs() {return coefs_;} //!< return coefs vector (write)
coefs() const130 const std::vector<complex_t>& coefs() const {return coefs_;} //!< return coefs vector (read)
leftPriority() const131 bool leftPriority() const {return leftPriority_;} //!< return the priority flag
leftPriority()132 bool& leftPriority() {return leftPriority_;} //!< return the priority flag
133 dimen_t diffOrder() const; //!< order of differential operator
134 bool normalRequired() const; //!< true if normal involved
135 bool xnormalRequired() const; //!< true if xnormal involved
136 bool ynormalRequired() const; //!< true if ynormal involved
137 bool xnormalRequiredByOpKernel() const; //!< true if normal is required by operator on Kernel
138 bool ynormalRequiredByOpKernel() const; //!< true if normal is required by operator on Kernel
extensionRequired() const139 bool extensionRequired() const //! true if operator involve non-tangential derivatives or space has no FE trace space
140 {return difOp_p->extensionRequired() || u_p->extensionRequired();}
141 bool elementRequired() const; //!< true if normal is required by operator
142
143
144 // utilities
145 number_t unknownDegree() const; //!< degree of unknown (basis function)
146 number_t degree() const; //!< degree of (derivative) unknown (basis function)
147 bool hasOperand() const; //!< true if operator involves left or right operand
148 bool hasFunction() const; //!< true if operator involves a user function or kernel
149 bool hasKernel() const; //!< true if operator involves a user Kernel
150 bool hasFun() const; //!< true if operator involves a user Function
isId() const151 bool isId() const //! true if id
152 {return difOp_p->type()==_id && leftOperand_p == 0 && rightOperand_p ==0;}
153 const Function* functionp() const; //!< direct access to Function if unique
154 const Kernel* kernelp() const; //!< direct access to Kernel
155 const OperatorOnKernel* opkernelp() const; //!< direct access to operator on Kernel
156 void changeKernel(Kernel*); //!< change Kernel in operator
setUnknown(const Unknown & u)157 void setUnknown(const Unknown& u) {u_p=&u;} //!< set (change) the unknowns
158
hasExtension() const159 bool hasExtension() const //! true if handles an extension
160 {if((leftOperand_p!=0 && leftOperand_p->hasExtension())
161 || (rightOperand_p!=0 && rightOperand_p->hasExtension())) return true;
162 return false;}
163
extension() const164 const Extension* extension() const //! return extension pointer, 0 if there is no extension
165 { const Extension* ext=0;
166 if(leftOperand_p!=0) ext=leftOperand_p->extension();
167 if(ext==0 && rightOperand_p!=0) ext=rightOperand_p->extension();
168 return ext;}
169
170 //transmit informations to user function or user kernel
setX(const Point & P) const171 void setX(const Point& P) const //! set X point if Function involved (usefull for Kernel)
172 {if(leftOperand_p!=0) leftOperand_p->setX(P);
173 if(rightOperand_p!=0) rightOperand_p->setX(P);}
setY(const Point & P) const174 void setY(const Point& P) const //! set Y point if Function involved (usefull for Kernel)
175 {if(leftOperand_p!=0) leftOperand_p->setY(P);
176 if(rightOperand_p!=0) rightOperand_p->setY(P);}
177
178 //Equation constructor (implemented in LcOperatorOnUnknown.cpp)
179 EssentialCondition operator = (const real_t &); //!< construct equation Op(u)=r
180 EssentialCondition operator = (const complex_t &); //!< construct equation Op(u)=c
181 EssentialCondition operator = (const Function &); //!< construct equation Op(u)=f
182 EssentialCondition operator = (const TermVector &); //!< construct equation Op(u)=TermVector (implemented in TermVector.cpp)
183
184 template <typename T,typename K>
185 void eval(const std::vector<K>&,
186 const std::vector< std::vector<K> >&,
187 dimen_t, Vector<T>&, dimen_t&, dimen_t&,
188 const Vector<real_t>* = 0 ) const; //!< evaluate operator from shape values (no function to evaluate)
189
190 // template <typename T,typename K>
191 // void eval(const Point&, const std::vector<K>&,
192 // const std::vector< std::vector<K> >&,
193 // dimen_t, Vector<T>&, dimen_t&, dimen_t&,
194 // const Vector<real_t>* = 0 ) const; //!< evaluate operator from point and shape values (function to evaluate)
195 template <typename T,typename K>
196 void eval(const Point& p, const std::vector<K>&,
197 const std::vector< std::vector<K> >&,
198 dimen_t, Vector<T>&, dimen_t&, dimen_t&,
199 const Vector<real_t>* =0, const ExtensionData* =0) const; //!< evaluate operator from point and shape values (function to evaluate)
200
201 template <typename T,typename K>
202 void eval(const Point&,const Point&, const std::vector<K>&,
203 const std::vector< std::vector<K> >&,
204 dimen_t, Vector<T>&, dimen_t&, dimen_t&,
205 const Vector<real_t>* = 0, const Vector<real_t>* = 0 ) const; //!< evaluate operator from point and shape values (kernel to evaluate)
206
207 void print(std::ostream&) const; //!< print attributes on stream
print(PrintStream & os) const208 void print(PrintStream& os) const {print(os.currentStream());}
209 void printsymbolic(std::ostream &) const; //!< print attributes in symbolic form on stream
printsymbolic(PrintStream & os) const210 void printsymbolic(PrintStream& os) const {printsymbolic(os.currentStream());}
211 string_t asString() const; //!< return as symbolic string
212 }; // end of class OperatorOnUnknown
213
214 // external functions
215 bool checkConsistancy(const OperatorOnUnknown&, AlgebraicOperator, const OperatorOnUnknown&); //!< check opu aop opv consistancy
216 bool operator==(const OperatorOnUnknown&, const OperatorOnUnknown&);//!< compare OperatorOnUnknown (same unknowns, same diff operators, same functions ...)
217 bool operator!=(const OperatorOnUnknown&, const OperatorOnUnknown&);//!< compare OperatorOnUnknown (same unknowns, same diff operators, same functions ...)
218
219 //=================================================================================
220 //@{
221 //! "differential" operators applied to unknown
222 OperatorOnUnknown& id(const Unknown&);
223 OperatorOnUnknown& d0(const Unknown&);
224 OperatorOnUnknown& dt(const Unknown&);
225 OperatorOnUnknown& d1(const Unknown&);
226 OperatorOnUnknown& dx(const Unknown&);
227 OperatorOnUnknown& d2(const Unknown&);
228 OperatorOnUnknown& dy(const Unknown&);
229 OperatorOnUnknown& d3(const Unknown&);
230 OperatorOnUnknown& dz(const Unknown&);
231 OperatorOnUnknown& grad(const Unknown&);
232 OperatorOnUnknown& nabla(const Unknown&);
233 OperatorOnUnknown& div(const Unknown&);
234 OperatorOnUnknown& curl(const Unknown&);
235 OperatorOnUnknown& rot(const Unknown&);
236 OperatorOnUnknown& gradS(const Unknown&);
237 OperatorOnUnknown& nablaS(const Unknown&);
238 OperatorOnUnknown& divS(const Unknown&);
239 OperatorOnUnknown& curlS(const Unknown&);
240 OperatorOnUnknown& rotS(const Unknown&);
241 OperatorOnUnknown& nx(const Unknown&);
242 OperatorOnUnknown& ndot(const Unknown&);
243 OperatorOnUnknown& ndotgrad(const Unknown&);
244 OperatorOnUnknown& ncrossgrad(const Unknown&);
245 OperatorOnUnknown& ncross(const Unknown&);
246 OperatorOnUnknown& ncrossncross(const Unknown&);
247 OperatorOnUnknown& ncrossntimes(const Unknown&);
248 OperatorOnUnknown& ncrossndot(const Unknown&);
249 OperatorOnUnknown& ntimesndot(const Unknown&);
250 OperatorOnUnknown& ndiv(const Unknown&);
251 OperatorOnUnknown& ncrosscurl(const Unknown&);
252 OperatorOnUnknown& epsilon(const Unknown&);
253 OperatorOnUnknown& epsilonR(const Unknown&);
254 OperatorOnUnknown& voigtToM(const Unknown&); //special operator : convert Voigt vector to tensor matrix
255 OperatorOnUnknown& gradG(const Unknown& un, const complex_t& ax, const complex_t& ay = complex_t(0., 0.),
256 const complex_t& az = complex_t(0., 0.) , const complex_t& at = complex_t(0., 0.));
257 OperatorOnUnknown& nablaG(const Unknown& un, const complex_t& ax, const complex_t& ay = complex_t(0., 0.),
258 const complex_t& az = complex_t(0., 0.) , const complex_t& at = complex_t(0., 0.));
259 OperatorOnUnknown& divG(const Unknown& un, const complex_t& ax, const complex_t& ay = complex_t(0., 0.),
260 const complex_t& az = complex_t(0., 0.) , const complex_t& at = complex_t(0., 0.));
261 OperatorOnUnknown& curlG(const Unknown& un, const complex_t& ax, const complex_t& ay = complex_t(0., 0.),
262 const complex_t& az = complex_t(0., 0.) , const complex_t& at = complex_t(0., 0.));
263 OperatorOnUnknown& rotG(const Unknown& un, const complex_t& ax, const complex_t& ay = complex_t(0., 0.),
264 const complex_t& az = complex_t(0., 0.) , const complex_t& at = complex_t(0., 0.));
265 OperatorOnUnknown& epsilonG(const Unknown& un, const complex_t& id, const complex_t& ax = complex_t(0., 0.),
266 const complex_t& ay = complex_t(0., 0.) , const complex_t& az = complex_t(0., 0.));
267 OperatorOnUnknown& mean(const Unknown&);
268 OperatorOnUnknown& jump(const Unknown&);
269 //@}
270
271 // for syntax difop&Unknown (for instance _grad&u means grad(u))
272 // OperatorOnUnknown& operator&(DiffOpType,const Unknown &);
273
274 //==========================================================================================
275 // equivalent differential operator with like mathematical forms (UnitaryVector and Unknown)
276 //==========================================================================================
277 OperatorOnUnknown& operator*(UnitaryVector, const Unknown&); //!< n*u same as nx(u)
278 OperatorOnUnknown& operator|(UnitaryVector, const Unknown&); //!< n|u same as ndot(u)
279 OperatorOnUnknown& operator^(UnitaryVector, const Unknown&); //!< n^u same as ncross(u)
280 OperatorOnUnknown& operator^(const Unknown&, UnitaryVector); //!< u^n same as -ncross(u)
281 OperatorOnUnknown& operator*(const Unknown&, UnitaryVector); //!< u*n same as nx(u)
282 OperatorOnUnknown& operator|(const Unknown&, UnitaryVector); //!< u|n same as ndot(u)
283 OperatorOnUnknown& operator|(UnitaryVector, OperatorOnUnknown&); //!< n|grad(u) same as ndotgrad(u)
284 OperatorOnUnknown& operator^(UnitaryVector, OperatorOnUnknown&); //!< n^(n^u) or n^curl(u) same as ncrossncross(u) or ncrosscurl(u)
285 OperatorOnUnknown& operator^(OperatorOnUnknown&, UnitaryVector); //!< -n^(n^u) or -n^curl(u) same as -ncrossncross(u) or -ncrosscurl(u)
286 OperatorOnUnknown& operator*(UnitaryVector, OperatorOnUnknown&); //!< n*div same as ndiv(u)
287 OperatorOnUnknown& operator|(OperatorOnUnknown&, UnitaryVector); //!< grad(u)|n same as ndotgrad(u)
288 OperatorOnUnknown& operator*(OperatorOnUnknown&, UnitaryVector); //!< div*n same as ndiv(u)
289
290 //==========================================================================================
291 // left and right operations with Function and Unknown
292 //==========================================================================================
293 OperatorOnUnknown& operator*(const Unknown&, const Function&); //!< u*F
294 OperatorOnUnknown& operator|(const Unknown&, const Function&); //!< u|F
295 OperatorOnUnknown& operator^(const Unknown&, const Function&); //!< u^F
296 OperatorOnUnknown& operator%(const Unknown&, const Function&); //!< u%F
297 OperatorOnUnknown& operator*(const Function&, const Unknown&); //!< F*u
298 OperatorOnUnknown& operator|(const Function&, const Unknown&); //!< F|u
299 OperatorOnUnknown& operator^(const Function&, const Unknown&); //!< F^u
300 OperatorOnUnknown& operator%(const Function&, const Unknown&); //!< F%u
301
302 //==========================================================================================
303 // left and right operations with OperatorOnFunction and Unknown
304 //==========================================================================================
305 OperatorOnUnknown& operator*(const Unknown&, const OperatorOnFunction&); //!< u*op(F)
306 OperatorOnUnknown& operator|(const Unknown&, const OperatorOnFunction&); //!< u|op(F)
307 OperatorOnUnknown& operator^(const Unknown&, const OperatorOnFunction&); //!< u^op(F)
308 OperatorOnUnknown& operator%(const Unknown&, const OperatorOnFunction&); //!< u%op(F)
309 OperatorOnUnknown& operator*(const OperatorOnFunction&, const Unknown&); //!< op(F)*u
310 OperatorOnUnknown& operator|(const OperatorOnFunction&, const Unknown&); //!< op(F)|u
311 OperatorOnUnknown& operator^(const OperatorOnFunction&, const Unknown&); //!< op(F)^u
312 OperatorOnUnknown& operator%(const OperatorOnFunction&, const Unknown&); //!< op(F)%u
313
314 //==========================================================================================
315 // left and right operations with C++ user functions and Unknown
316 // user shortcuts, avoid using explicit Function encapsulation
317 // only function, no kernel (see below)
318 //==========================================================================================
319 template <typename T>
operator *(const Unknown & un,T (fun)(const Point &,Parameters &))320 OperatorOnUnknown& operator*(const Unknown& un, T(fun)(const Point&, Parameters&)) //! u * function(Point,Parameters)
321 {
322 Function& f = *(new Function(*fun)); // create a Function object on heap memory
323 return un * f;
324 }
325 template <typename T>
operator *(T (fun)(const Point &,Parameters &),const Unknown & un)326 OperatorOnUnknown& operator*(T(fun)(const Point&, Parameters&), const Unknown& un) //! function(Point,Parameters) * u
327 {
328 Function& f = *(new Function(*fun)); // create a Function object on heap memory
329 return f * un;
330 }
331 template <typename T>
operator |(const Unknown & un,T (fun)(const Point &,Parameters &))332 OperatorOnUnknown& operator|(const Unknown& un, T(fun)(const Point&, Parameters&)) //! u | function(Point,Parameters)
333 {
334 Function& f = *(new Function(*fun)); // create a Function object on heap memory
335 return un | f;
336 }
337 template <typename T>
operator |(T (fun)(const Point &,Parameters &),const Unknown & un)338 OperatorOnUnknown& operator|(T(fun)(const Point&, Parameters&), const Unknown& un) //! function(Point,Parameters) | u
339 {
340 Function& f = *(new Function(*fun)); // create a Function object on heap memory
341 return f | un;
342 }
343 template <typename T>
operator ^(const Unknown & un,T (fun)(const Point &,Parameters &))344 OperatorOnUnknown& operator^(const Unknown& un, T(fun)(const Point&, Parameters&)) //! u ^ function(Point,Parameters)
345 {
346 Function& f = *(new Function(*fun)); // create a Function object on heap memory
347 return un ^ f;
348 }
349 template <typename T>
operator ^(T (fun)(const Point &,Parameters &),const Unknown & un)350 OperatorOnUnknown& operator^(T(fun)(const Point&, Parameters&), const Unknown& un) //! function(Point,Parameters) ^ u
351 {
352 Function& f = *(new Function(*fun)); // create a Function object on heap memory
353 return f ^ un;
354 }
355 template <typename T>
operator %(const Unknown & un,T (fun)(const Point &,Parameters &))356 OperatorOnUnknown& operator%(const Unknown& un, T(fun)(const Point&, Parameters&)) //! u % function(Point,Parameters)
357 {
358 Function& f = *(new Function(*fun)); // create a Function object on heap memory
359 return un % f;
360 }
361 template <typename T>
operator %(T (fun)(const Point &,Parameters &),const Unknown & un)362 OperatorOnUnknown& operator%(T(fun)(const Point&, Parameters&), const Unknown& un) //! function(Point,Parameters) % u
363 {
364 Function& f = *(new Function(*fun)); // create a Function object on heap memory
365 return f % un;
366 }
367
368 // same for vector function (Vector<Point> as argument)
369 template <typename T>
operator *(const Unknown & un,T (fun)(const Vector<Point> &,Parameters &))370 OperatorOnUnknown& operator*(const Unknown& un, T(fun)(const Vector<Point>&, Parameters&)) //! u * function(Vector<Point>,Parameters)
371 {
372 Function& f = *(new Function(*fun)); // create a Function object on heap memory
373 return un * f;
374 }
375 template <typename T>
operator *(T (fun)(const Vector<Point> &,Parameters &),const Unknown & un)376 OperatorOnUnknown& operator*(T(fun)(const Vector<Point>&, Parameters&), const Unknown& un) //! function(Vector<Point>,Parameters) * u
377 {
378 Function& f = *(new Function(*fun)); // create a Function object on heap memory
379 return f * un;
380 }
381 template <typename T>
operator |(const Unknown & un,T (fun)(const Vector<Point> &,Parameters &))382 OperatorOnUnknown& operator|(const Unknown& un, T(fun)(const Vector<Point>&, Parameters&)) //! u | function(Vector<Point>,Parameters)
383 {
384 Function& f = *(new Function(*fun)); // create a Function object on heap memory
385 return un | f;
386 }
387 template <typename T>
operator |(T (fun)(const Vector<Point> &,Parameters &),const Unknown & un)388 OperatorOnUnknown& operator|(T(fun)(const Vector<Point>&, Parameters&), const Unknown& un) //! function(Vector<Point>,Parameters) | u
389 {
390 Function& f = *(new Function(*fun)); // create a Function object on heap memory
391 return f | un;
392 }
393 template <typename T>
operator ^(const Unknown & un,T (fun)(const Vector<Point> &,Parameters &))394 OperatorOnUnknown& operator^(const Unknown& un, T(fun)(const Vector<Point>&, Parameters&)) //! u ^ function(Vector<Point>,Parameters)
395 {
396 Function& f = *(new Function(*fun)); // create a Function object on heap memory
397 return un ^ f;
398 }
399 template <typename T>
operator ^(T (fun)(const Vector<Point> &,Parameters &),const Unknown & un)400 OperatorOnUnknown& operator^(T(fun)(const Vector<Point>&, Parameters&), const Unknown& un) //! function(Vector<Point>,Parameters) ^ u
401 {
402 Function& f = *(new Function(*fun)); // create a Function object on heap memory
403 return f ^ un;
404 }
405 template <typename T>
operator %(const Unknown & un,T (fun)(const Vector<Point> &,Parameters &))406 OperatorOnUnknown& operator%(const Unknown& un, T(fun)(const Vector<Point>&, Parameters&)) //! u % function(Vector<Point>,Parameters)
407 {
408 Function& f = *(new Function(*fun)); // create a Function object on heap memory
409 return un % f;
410 }
411 template <typename T>
operator %(T (fun)(const Vector<Point> &,Parameters &),const Unknown & un)412 OperatorOnUnknown& operator%(T(fun)(const Vector<Point>&, Parameters&), const Unknown& un) //! function(Vector<Point>,Parameters) % u
413 {
414 Function& f = *(new Function(*fun)); // create a Function object on heap memory
415 return f % un;
416 }
417
418 //==========================================================================================
419 // left and right operations with Value and Unknown
420 //==========================================================================================
421 OperatorOnUnknown& operator*(const Unknown&, const Value&); //!< u*val
422 OperatorOnUnknown& operator|(const Unknown&, const Value&); //!< u|val
423 OperatorOnUnknown& operator^(const Unknown&, const Value&); //!< u^val
424 OperatorOnUnknown& operator%(const Unknown&, const Value&); //!< u%val
425 OperatorOnUnknown& operator*(const Value&, const Unknown&); //!< val*u
426 OperatorOnUnknown& operator|(const Value&, const Unknown&); //!< val|u
427 OperatorOnUnknown& operator^(const Value&, const Unknown&); //!< val^u
428 OperatorOnUnknown& operator%(const Value&, const Unknown&); //!< val%u
429
430 //==========================================================================================
431 // left and right operations with any value and Unknown
432 // to avoid intricated conflict we do not use templated operators
433 //==========================================================================================
434 template<typename T>
fromUnknownVal(const Unknown & un,const T & val,AlgebraicOperator aop)435 OperatorOnUnknown& fromUnknownVal(const Unknown& un, const T& val, AlgebraicOperator aop)
436 {
437 OperatorOnUnknown* opu = new OperatorOnUnknown(&un, _id);
438 updateRight(*opu, val, aop);
439 return *opu;
440 }
441
442 template<typename T>
fromValUnknown(const Unknown & un,const T & val,AlgebraicOperator aop)443 OperatorOnUnknown& fromValUnknown(const Unknown& un, const T& val, AlgebraicOperator aop)
444 {
445 OperatorOnUnknown* opu = new OperatorOnUnknown(&un, _id);
446 updateLeft(*opu, val, aop);
447 return *opu;
448 }
449
450 //unknown op r or r op unknown
451 OperatorOnUnknown& operator*(const Unknown&, const real_t&); //!< u*r
452 OperatorOnUnknown& operator|(const Unknown&, const real_t&); //!< u|r
453 OperatorOnUnknown& operator^(const Unknown&, const real_t&); //!< u^r
454 OperatorOnUnknown& operator%(const Unknown&, const real_t&); //!< u%r
455 OperatorOnUnknown& operator*(const real_t&, const Unknown&); //!< r*u
456 OperatorOnUnknown& operator|(const real_t&, const Unknown&); //!< r|u
457 OperatorOnUnknown& operator^(const real_t&, const Unknown&); //!< r^u
458 OperatorOnUnknown& operator%(const real_t&, const Unknown&); //!< r%u
459
460 //unknown op c or c op unknown
461 OperatorOnUnknown& operator*(const Unknown&, const complex_t&); //!< u*r
462 OperatorOnUnknown& operator|(const Unknown&, const complex_t&); //!< u|c
463 OperatorOnUnknown& operator^(const Unknown&, const complex_t&); //!< u^c
464 OperatorOnUnknown& operator%(const Unknown&, const complex_t&); //!< u%c
465 OperatorOnUnknown& operator*(const complex_t&, const Unknown&); //!< c*u
466 OperatorOnUnknown& operator|(const complex_t&, const Unknown&); //!< c|u
467 OperatorOnUnknown& operator^(const complex_t&, const Unknown&); //!< c^u
468 OperatorOnUnknown& operator%(const complex_t&, const Unknown&); //!< c%u
469
470 //unknown op vector or vector op unknown
471 template<typename T>
operator *(const Unknown & un,const Vector<T> & val)472 OperatorOnUnknown& operator*(const Unknown& un, const Vector<T>& val) //! u*Vector
473 {return fromUnknownVal(un,val,_product);}
474 template<typename T>
operator |(const Unknown & un,const Vector<T> & val)475 OperatorOnUnknown& operator|(const Unknown& un, const Vector<T>& val) //! u|Vector
476 {return fromUnknownVal(un,val,_innerProduct);}
477 template<typename T>
operator ^(const Unknown & un,const Vector<T> & val)478 OperatorOnUnknown& operator^(const Unknown& un, const Vector<T>& val) //! u^Vector
479 {return fromUnknownVal(un,val,_crossProduct);}
480 template<typename T>
operator %(const Unknown & un,const Vector<T> & val)481 OperatorOnUnknown& operator%(const Unknown& un, const Vector<T>& val) //! u%Vector
482 {return fromUnknownVal(un,val,_contractedProduct);}
483 template<typename T>
operator *(const Vector<T> & val,const Unknown & un)484 OperatorOnUnknown& operator*(const Vector<T>& val, const Unknown& un) //! Vector*u
485 {return fromValUnknown(un,val,_product);}
486 template<typename T>
operator |(const Vector<T> & val,const Unknown & un)487 OperatorOnUnknown& operator|(const Vector<T>& val, const Unknown& un) //! Vector|u
488 {return fromValUnknown(un,val,_innerProduct);}
489 template<typename T>
operator ^(const Vector<T> & val,const Unknown & un)490 OperatorOnUnknown& operator^(const Vector<T>& val, const Unknown& un) //! Vector^u
491 {return fromValUnknown(un,val,_crossProduct);}
492 template<typename T>
operator %(const Vector<T> & val,const Unknown & un)493 OperatorOnUnknown& operator%(const Vector<T>& val, const Unknown& un) //! Vector%u
494 {return fromValUnknown(un,val,_contractedProduct);}
495
496 //unknown op matrix or matrix op unknown
497 template<typename T>
operator *(const Unknown & un,const Matrix<T> & val)498 OperatorOnUnknown& operator*(const Unknown& un, const Matrix<T>& val) //! u*Matrix
499 {return fromUnknownVal(un,val,_product);}
500 template<typename T>
operator |(const Unknown & un,const Matrix<T> & val)501 OperatorOnUnknown& operator|(const Unknown& un, const Matrix<T>& val) //! u|Matrix
502 {return fromUnknownVal(un,val,_innerProduct);}
503 template<typename T>
operator ^(const Unknown & un,const Matrix<T> & val)504 OperatorOnUnknown& operator^(const Unknown& un, const Matrix<T>& val) //! u^Matrix
505 {return fromUnknownVal(un,val,_crossProduct);}
506 template<typename T>
operator %(const Unknown & un,const Matrix<T> & val)507 OperatorOnUnknown& operator%(const Unknown& un, const Matrix<T>& val) //! u%Matrix
508 {return fromUnknownVal(un,val,_contractedProduct);}
509 template<typename T>
operator *(const Matrix<T> & val,const Unknown & un)510 OperatorOnUnknown& operator*(const Matrix<T>& val, const Unknown& un) //! Matrix*u
511 {return fromValUnknown(un,val,_product);}
512 template<typename T>
operator |(const Matrix<T> & val,const Unknown & un)513 OperatorOnUnknown& operator|(const Matrix<T>& val, const Unknown& un) //! Matrix|u
514 {return fromValUnknown(un,val,_innerProduct);}
515 template<typename T>
operator ^(const Matrix<T> & val,const Unknown & un)516 OperatorOnUnknown& operator^(const Matrix<T>& val, const Unknown& un) //! Matrix^u
517 {return fromValUnknown(un,val,_crossProduct);}
518 template<typename T>
operator %(const Matrix<T> & val,const Unknown & un)519 OperatorOnUnknown& operator%(const Matrix<T>& val, const Unknown& un) //! Matrix%u
520 {return fromValUnknown(un,val,_contractedProduct);}
521
522 //==========================================================================================
523 // left and right operation with Function : op(u) aop Function or Function aop op(u)
524 //==========================================================================================
525 OperatorOnUnknown& updateRight(OperatorOnUnknown&, const Function&, AlgebraicOperator); //!< update Op(u)*F operation
526 OperatorOnUnknown& updateLeft(OperatorOnUnknown&, const Function&, AlgebraicOperator); //!< update F*Op(u) operation
527 OperatorOnUnknown& operator*(OperatorOnUnknown&, const Function&); //!< product syntax Op(u)*F
528 OperatorOnUnknown& operator|(OperatorOnUnknown&, const Function&); //!< innerproduct syntax Op(u)|F
529 OperatorOnUnknown& operator^(OperatorOnUnknown&, const Function&); //!< cross product syntax Op(u)^F
530 OperatorOnUnknown& operator%(OperatorOnUnknown&, const Function&); //!< contracted product syntax Op(u)%F
531 OperatorOnUnknown& operator*(const Function&, OperatorOnUnknown&); //!< product syntax F*Op(u)
532 OperatorOnUnknown& operator|(const Function&, OperatorOnUnknown&); //!< innerproduct syntax F|Op(u)
533 OperatorOnUnknown& operator^(const Function&, OperatorOnUnknown&); //!< cross product syntax F^Op(u)
534 OperatorOnUnknown& operator%(const Function&, OperatorOnUnknown&); //!< contracted product syntax F%Op(u)
535
536 //==========================================================================================
537 // left and right operation with Function : op(u) aop op(F) or op(F) aop op(u)
538 //==========================================================================================
539 OperatorOnUnknown& updateRight(OperatorOnUnknown&, const OperatorOnFunction&, AlgebraicOperator); //!< update Op(u)*op(F) operation
540 OperatorOnUnknown& updateLeft(OperatorOnUnknown&, const OperatorOnFunction&, AlgebraicOperator); //!< update op(F)*Op(u) operation
541 OperatorOnUnknown& operator*(OperatorOnUnknown&, const OperatorOnFunction&); //!< product syntax Op(u)*op(F)
542 OperatorOnUnknown& operator|(OperatorOnUnknown&, const OperatorOnFunction&); //!< innerproduct syntax Op(u)|op(F)
543 OperatorOnUnknown& operator^(OperatorOnUnknown&, const OperatorOnFunction&); //!< cross product syntax Op(u)^op(F)
544 OperatorOnUnknown& operator%(OperatorOnUnknown&, const OperatorOnFunction&); //!< contracted product syntax Op(u)%op(F)
545 OperatorOnUnknown& operator*(const OperatorOnFunction&, OperatorOnUnknown&); //!< product syntax op(F)*Op(u)
546 OperatorOnUnknown& operator|(const OperatorOnFunction&, OperatorOnUnknown&); //!< innerproduct syntax op(F)|Op(u)
547 OperatorOnUnknown& operator^(const OperatorOnFunction&, OperatorOnUnknown&); //!< cross product syntax op(F)^Op(u)
548 OperatorOnUnknown& operator%(const OperatorOnFunction&, OperatorOnUnknown&); //!< contracted product syntax op(F)F%Op(u)
549
550 //==========================================================================================
551 // left and right operation with Value : op(u) aop Value or Value aop op(u)
552 //==========================================================================================
553 OperatorOnUnknown& updateRight(OperatorOnUnknown&, const Value&, AlgebraicOperator); //!< update Op(u)*V operation
554 OperatorOnUnknown& updateLeft(OperatorOnUnknown&, const Value&, AlgebraicOperator); //!< update V*Op(u) operation
555 OperatorOnUnknown& operator*(OperatorOnUnknown&, const Value&); //!< product syntax Op(u)*V
556 OperatorOnUnknown& operator|(OperatorOnUnknown&, const Value&); //!< innerproduct syntax Op(u)|V
557 OperatorOnUnknown& operator^(OperatorOnUnknown&, const Value&); //!< cross product syntax Op(u)^V
558 OperatorOnUnknown& operator%(OperatorOnUnknown&, const Value&); //!< contracted product syntax Op(u)%V
559 OperatorOnUnknown& operator*(const Value&, OperatorOnUnknown&); //!< product syntax V*Op(u)
560 OperatorOnUnknown& operator|(const Value&, OperatorOnUnknown&); //!< innerproduct syntax V|Op(u)
561 OperatorOnUnknown& operator^(const Value&, OperatorOnUnknown&); //!< cross product syntax V^Op(u)
562 OperatorOnUnknown& operator%(const Value&, OperatorOnUnknown&); //!< contracted product syntax V%Op(u)
563
564 //==========================================================================================
565 // update left and right operation with anything : op(u) aop anything or anything aop op(u)
566 //==========================================================================================
updateRight(OperatorOnUnknown & opu,const T & val,AlgebraicOperator o)567 template<typename T> OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const T& val, AlgebraicOperator o)
568 {
569 Value::checkTypeInList(val); // check the type of val
570 if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
571 opu.rightOperand() = new Operand(val, o);
572 // opu.updateReturnedType(o, opu.rightOperand()->value().valueType(),
573 // opu.rightOperand()->value().strucType(), opu.rightOperand()->value().dims(), false);
574 if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
575 opu.setStructure();
576 return opu;
577 }
578
updateLeft(OperatorOnUnknown & opu,const T & val,AlgebraicOperator o)579 template<typename T> OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const T& val, AlgebraicOperator o)
580 {
581 Value::checkTypeInList(val); // check the type of val
582 if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
583 opu.leftOperand() = new Operand(val, o);
584 // opu.updateReturnedType(o, opu.leftOperand()->value().valueType(),
585 // opu.leftOperand()->value().strucType(), opu.leftOperand()->value().dims(),true);
586 if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
587 opu.setStructure();
588 return opu;
589 }
590
591 // left and right operation with real scalar : op(u) aop r or r aop op(u)
592 OperatorOnUnknown& operator*(OperatorOnUnknown&, const real_t&);
593 OperatorOnUnknown& operator|(OperatorOnUnknown&, const real_t&);
594 OperatorOnUnknown& operator^(OperatorOnUnknown&, const real_t&);
595 OperatorOnUnknown& operator%(OperatorOnUnknown&, const real_t&);
596 OperatorOnUnknown& operator*(const real_t&, OperatorOnUnknown&);
597 OperatorOnUnknown& operator|(const real_t&, OperatorOnUnknown&);
598 OperatorOnUnknown& operator^(const real_t&, OperatorOnUnknown&);
599 OperatorOnUnknown& operator%(const real_t&, OperatorOnUnknown&);
600
601 // left and right operation with complex scalar : op(u) aop c or c aop op(u)
602 OperatorOnUnknown& operator*(OperatorOnUnknown&, const complex_t&);
603 OperatorOnUnknown& operator|(OperatorOnUnknown&, const complex_t&);
604 OperatorOnUnknown& operator^(OperatorOnUnknown&, const complex_t&);
605 OperatorOnUnknown& operator%(OperatorOnUnknown&, const complex_t&);
606 OperatorOnUnknown& operator*(const complex_t&, OperatorOnUnknown&);
607 OperatorOnUnknown& operator|(const complex_t&, OperatorOnUnknown&);
608 OperatorOnUnknown& operator^(const complex_t&, OperatorOnUnknown&);
609 OperatorOnUnknown& operator%(const complex_t&, OperatorOnUnknown&);
610
611 // left and right operation with vector<T> : op(u) aop vector or vector aop op(u)
operator *(OperatorOnUnknown & opu,const Vector<T> & val)612 template<typename T> OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const Vector<T>& val)
613 {return updateRight(opu, val, _product);}
operator |(OperatorOnUnknown & opu,const Vector<T> & val)614 template<typename T> OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const Vector<T>& val)
615 {return updateRight(opu, val, _innerProduct);}
operator ^(OperatorOnUnknown & opu,const Vector<T> & val)616 template<typename T> OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const Vector<T>& val)
617 {return updateRight(opu, val, _crossProduct);}
operator %(OperatorOnUnknown & opu,const Vector<T> & val)618 template<typename T> OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const Vector<T>& val)
619 {return updateRight(opu, val, _contractedProduct);}
operator *(const Vector<T> & val,OperatorOnUnknown & opu)620 template<typename T> OperatorOnUnknown& operator*(const Vector<T>& val, OperatorOnUnknown& opu)
621 {return updateLeft(opu, val, _product);}
operator |(const Vector<T> & val,OperatorOnUnknown & opu)622 template<typename T> OperatorOnUnknown& operator|(const Vector<T>& val, OperatorOnUnknown& opu)
623 {return updateLeft(opu, val, _innerProduct);}
operator ^(const Vector<T> & val,OperatorOnUnknown & opu)624 template<typename T> OperatorOnUnknown& operator^(const Vector<T>& val, OperatorOnUnknown& opu)
625 {return updateLeft(opu, val, _crossProduct);}
operator %(const Vector<T> & val,OperatorOnUnknown & opu)626 template<typename T> OperatorOnUnknown& operator%(const Vector<T>& val, OperatorOnUnknown& opu)
627 {return updateLeft(opu, val, _contractedProduct);}
628
629 // left and right operation with Matrix<T> : op(u) aop matrix or matrix aop op(u)
operator *(OperatorOnUnknown & opu,const Matrix<T> & val)630 template<typename T> OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const Matrix<T>& val)
631 {return updateRight(opu, val, _product);}
operator |(OperatorOnUnknown & opu,const Matrix<T> & val)632 template<typename T> OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const Matrix<T>& val)
633 {return updateRight(opu, val, _innerProduct);}
operator ^(OperatorOnUnknown & opu,const Matrix<T> & val)634 template<typename T> OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const Matrix<T>& val)
635 {return updateRight(opu, val, _crossProduct);}
operator %(OperatorOnUnknown & opu,const Matrix<T> & val)636 template<typename T> OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const Matrix<T>& val)
637 {return updateRight(opu, val, _contractedProduct);}
operator *(const Matrix<T> & val,OperatorOnUnknown & opu)638 template<typename T> OperatorOnUnknown& operator*(const Matrix<T>& val, OperatorOnUnknown& opu)
639 {return updateLeft(opu, val, _product);}
operator |(const Matrix<T> & val,OperatorOnUnknown & opu)640 template<typename T> OperatorOnUnknown& operator|(const Matrix<T>& val, OperatorOnUnknown& opu)
641 {return updateLeft(opu, val, _innerProduct);}
operator ^(const Matrix<T> & val,OperatorOnUnknown & opu)642 template<typename T> OperatorOnUnknown& operator^(const Matrix<T>& val, OperatorOnUnknown& opu)
643 {return updateLeft(opu, val, _crossProduct);}
operator %(const Matrix<T> & val,OperatorOnUnknown & opu)644 template<typename T> OperatorOnUnknown& operator%(const Matrix<T>& val, OperatorOnUnknown& opu)
645 {return updateLeft(opu, val, _contractedProduct);}
646
647 // left and right operation with anything: op(u) aop anything or anything aop op(u)
648 //template<typename T> OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const T& val)
649 //{return updateRight(opu, val, _product);}
650 //template<typename T> OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const T& val)
651 //{return updateRight(opu, val, _innerProduct);}
652 //template<typename T> OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const T& val)
653 //{return updateRight(opu, val, _crossProduct);}
654 //template<typename T> OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const T& val)
655 //{return updateRight(opu, val, _contractedProduct);}
656 //template<typename T> OperatorOnUnknown& operator*(const T& val, OperatorOnUnknown& opu)
657 //{return updateLeft(opu, val, _product);}
658 //template<typename T> OperatorOnUnknown& operator|(const T& val, OperatorOnUnknown& opu)
659 //{return updateLeft(opu, val, _innerProduct);}
660 //template<typename T> OperatorOnUnknown& operator^(const T& val, OperatorOnUnknown& opu)
661 //{return updateLeft(opu, val, _crossProduct);}
662 //template<typename T> OperatorOnUnknown& operator%(const T& val, OperatorOnUnknown& opu)
663 //{return updateLeft(opu, val, _contractedProduct);}
664
665 //==========================================================================================
666 // left and right operations with C++ user functions and OperatorOnUnknown
667 // user shortcuts, avoid using explicit Function encapsulation
668 // only function, no kernel (see KernelOperatorOnUnknowns)
669 //==========================================================================================
670 template <typename T>
operator *(OperatorOnUnknown & opu,T (fun)(const Point &,Parameters &))671 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, T(fun)(const Point&, Parameters&)) //! opu * function(Point,Parameters)
672 {
673 Function& f = *(new Function(*fun)); // create a Function object on heap memory
674 return opu * f;
675 }
676 template <typename T>
operator *(T (fun)(const Point &,Parameters &),OperatorOnUnknown & opu)677 OperatorOnUnknown& operator*(T(fun)(const Point&, Parameters&), OperatorOnUnknown& opu) //! function(Point,Parameters) * opu
678 {
679 Function& f = *(new Function(*fun)); // create a Function object on heap memory
680 return f * opu;
681 }
682 template <typename T>
operator |(OperatorOnUnknown & opu,T (fun)(const Point &,Parameters &))683 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, T(fun)(const Point&, Parameters&)) //! opu | function(Point,Parameters)
684 {
685 Function& f = *(new Function(*fun)); // create a Function object on heap memory
686 return opu | f;
687 }
688 template <typename T>
operator |(T (fun)(const Point &,Parameters &),OperatorOnUnknown & opu)689 OperatorOnUnknown& operator|(T(fun)(const Point&, Parameters&), OperatorOnUnknown& opu) //! function(Point,Parameters) | opu
690 {
691 Function& f = *(new Function(*fun)); // create a Function object on heap memory
692 return f | opu;
693 }
694 template <typename T>
operator ^(OperatorOnUnknown & opu,T (fun)(const Point &,Parameters &))695 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, T(fun)(const Point&, Parameters&)) //! opu ^ function(Point,Parameters)
696 {
697 Function& f = *(new Function(*fun)); // create a Function object on heap memory
698 return opu ^ f;
699 }
700 template <typename T>
operator ^(T (fun)(const Point &,Parameters &),OperatorOnUnknown & opu)701 OperatorOnUnknown& operator^(T(fun)(const Point&, Parameters&), OperatorOnUnknown& opu) //! function(Point,Parameters) ^ opu
702 {
703 Function& f = *(new Function(*fun)); // create a Function object on heap memory
704 return f ^ opu;
705 }
706 template <typename T>
operator %(OperatorOnUnknown & opu,T (fun)(const Point &,Parameters &))707 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, T(fun)(const Point&, Parameters&)) //! opu % function(Point,Parameters)
708 {
709 Function& f = *(new Function(*fun)); // create a Function object on heap memory
710 return opu % f;
711 }
712 template <typename T>
operator %(T (fun)(const Point &,Parameters &),OperatorOnUnknown & opu)713 OperatorOnUnknown& operator%(T(fun)(const Point&, Parameters&), OperatorOnUnknown& opu) //! function(Point,Parameters) % opu
714 {
715 Function& f = *(new Function(*fun)); // create a Function object on heap memory
716 return f % opu;
717 }
718 template <typename T>
operator *(OperatorOnUnknown & opu,T (fun)(const Vector<Point> &,Parameters &))719 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, T(fun)(const Vector<Point>&, Parameters&)) //! opu * function(Vector<Point>,Parameters)
720 {
721 Function& f = *(new Function(*fun)); // create a Function object on heap memory
722 return opu * f;
723 }
724 template <typename T>
operator *(T (fun)(const Vector<Point> &,Parameters &),OperatorOnUnknown & opu)725 OperatorOnUnknown& operator*(T(fun)(const Vector<Point>&, Parameters&), OperatorOnUnknown& opu) //! function(Vector<Point>,Parameters) * opu
726 {
727 Function& f = *(new Function(*fun)); // create a Function object on heap memory
728 return f * opu;
729 }
730 template <typename T>
operator |(OperatorOnUnknown & opu,T (fun)(const Vector<Point> &,Parameters &))731 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, T(fun)(const Vector<Point>&, Parameters&)) //! opu | function(Vector<Point>,Parameters)
732 {
733 Function& f = *(new Function(*fun)); // create a Function object on heap memory
734 return opu | f;
735 }
736 template <typename T>
operator |(T (fun)(const Vector<Point> &,Parameters &),OperatorOnUnknown & opu)737 OperatorOnUnknown& operator|(T(fun)(const Vector<Point>&, Parameters&), OperatorOnUnknown& opu) //! function(Vector<Point>,Parameters) | opu
738 {
739 Function& f = *(new Function(*fun)); // create a Function object on heap memory
740 return f | opu;
741 }
742 template <typename T>
operator ^(OperatorOnUnknown & opu,T (fun)(const Vector<Point> &,Parameters &))743 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, T(fun)(const Vector<Point>&, Parameters&)) //! opu ^ function(Vector<Point>,Parameters)
744 {
745 Function& f = *(new Function(*fun)); // create a Function object on heap memory
746 return opu ^ f;
747 }
748 template <typename T>
operator ^(T (fun)(const Vector<Point> &,Parameters &),OperatorOnUnknown & opu)749 OperatorOnUnknown& operator^(T(fun)(const Vector<Point>&, Parameters&), OperatorOnUnknown& opu) //! function(Vector<Point>,Parameters) ^ opu
750 {
751 Function& f = *(new Function(*fun)); // create a Function object on heap memory
752 return f ^ opu;
753 }
754 template <typename T>
operator %(OperatorOnUnknown & opu,T (fun)(const Vector<Point> &,Parameters &))755 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, T(fun)(const Vector<Point>&, Parameters&)) //! opu % function(Vector<Point>,Parameters)
756 {
757 Function& f = *(new Function(*fun)); // create a Function object on heap memory
758 return opu % f;
759 }
760 template <typename T>
operator %(T (fun)(const Vector<Point> &,Parameters &),OperatorOnUnknown & opu)761 OperatorOnUnknown& operator%(T(fun)(const Vector<Point>&, Parameters&), OperatorOnUnknown& opu) //! function(Vector<Point>,Parameters) % opu
762 {
763 Function& f = *(new Function(*fun)); // create a Function object on heap memory
764 return f % opu;
765 }
766
767 //------------------------------------------------------------------------------------------------------------
768 /*!
769 evaluate operator from shape values and its derivatives
770 sv : shape values
771 dsv : shape derivatives
772 dimFun : dimension of shape function
773 val : operator values as a vector stored as following:
774 O1 O2 .... Os s number of shape functions
775 Ok may be
776 - a scalar
777 - a vector O_k=[Ok1, Ok2, ..., Okn]
778 - a matrix O_k=[Ok11, Ok12, ..., Ok1n, Ok21, Ok22, ..., Ok2n, ... Okm1, Okm2, ..., Okmn]
779 d : block size
780 m : number of matrix rows when returning some matrices
781 np : pointer to the normal vector, 0 if no available normal vector
782
783 T type of result val has to be consistent (no check)
784 K type of shape functions
785
786 There are 3 versions of eval corresponding to the cases:
787 - no function to evaluate
788 - function to evaluate at a point
789 - kernel to evaluate at a pair of points
790 */
791 //------------------------------------------------------------------------------------------------------------
792 template <typename T,typename K>
eval(const std::vector<K> & sv,const std::vector<std::vector<K>> & dsv,dimen_t dimFun,Vector<T> & val,dimen_t & d,dimen_t & m,const Vector<real_t> * np) const793 void OperatorOnUnknown::eval(const std::vector<K>& sv,
794 const std::vector< std::vector<K> >& dsv,
795 dimen_t dimFun, Vector<T>& val, dimen_t& d, dimen_t& m,
796 const Vector<real_t>* np) const
797 {
798 // compute differential operator
799 number_t nbw=sv.size()/dimFun;
800 d=dimFun;
801 Vector<K> difop;
802 difOp_p->eval(sv, dsv, d, m, difop, np, coefs_);
803
804 if(leftOperand_p == 0 && rightOperand_p == 0) {val=difop; return;}
805 if(leftOperand_p != 0 && rightOperand_p == 0) {val=leftOperand_p->leftEval<T>(difop, d, m, nbw);return;}
806 if(leftOperand_p == 0 && rightOperand_p != 0) {val=rightOperand_p->rightEval<T>(difop, d, m, nbw);return;}
807 val=leftOperand_p->leftEval<T>(rightOperand_p->rightEval<T>(difop, d, m, nbw), d, m, nbw);
808 }
809
810 // // evaluate operator from a point p and shape values sv and its derivatives dsv
811 // // (dimFun is the dimension of shape functions)
812 // // template <typename T,typename K>
813 // void OperatorOnUnknown::eval(const Point& p, const std::vector<K>& sv,
814 // const std::vector< std::vector<K> >& dsv,
815 // dimen_t dimFun, Vector<T>& val, dimen_t& d, dimen_t& m,
816 // const Vector<real_t>* np) const
817 // {
818 // // compute differential operator
819 // number_t nbw=sv.size()/dimFun;
820 // d=dimFun;
821 // Vector<K> difop;
822 // difOp_p->eval(sv, dsv, d, m, difop, np, coefs_);
823
824 // //compute operator on unknown
825 // if(leftOperand_p == 0 && rightOperand_p == 0) {val=difop; return;}
826 // if(leftOperand_p != 0 && rightOperand_p == 0) {val=leftOperand_p->leftEval<T>(p, difop, d, m, nbw, np); return;}
827 // if(leftOperand_p == 0 && rightOperand_p != 0) {val=rightOperand_p->rightEval<T>(p, difop, d, m, nbw, np); return;}
828 // if(leftOperand_p->isFunction()) val=leftOperand_p->leftEval<T>(p, difop, d, m, nbw, np);
829 // else val=leftOperand_p->leftEval<T>(difop,d,m,nbw);
830 // if(rightOperand_p->isFunction()) val=rightOperand_p->rightEval<T>(p, val, d, m, nbw, np);
831 // else val=rightOperand_p->rightEval<T>(val, d, m, nbw);
832 // }
833
834 //evaluate operator from a point p and shape values sv and its derivatives dsv when extension involved
835 // ExtensionData embends some useful data related to the extension (shapevalues, nodes on current element)
836 //(dimFun is the dimension of shape functions)
837 template <typename T,typename K>
eval(const Point & p,const std::vector<K> & sv,const std::vector<std::vector<K>> & dsv,dimen_t dimFun,Vector<T> & val,dimen_t & d,dimen_t & m,const Vector<real_t> * np,const ExtensionData * extdata) const838 void OperatorOnUnknown::eval(const Point& p, const std::vector<K>& sv,
839 const std::vector< std::vector<K> >& dsv,
840 dimen_t dimFun, Vector<T>& val, dimen_t& d, dimen_t& m,
841 const Vector<real_t>* np, const ExtensionData* extdata) const
842 {
843 // compute differential operator
844 number_t nbw=sv.size()/dimFun;
845 d=dimFun;
846 Vector<K> difop;
847 difOp_p->eval(sv, dsv, d, m, difop, np, coefs_);
848
849 //compute operator on unknown
850 if(leftOperand_p == 0)
851 {
852 if(rightOperand_p == 0) {val=difop; return;}
853 val=rightOperand_p->rightEval<T>(p, difop, d, m, nbw, np, extdata);
854 return;
855 }
856 else //leftoperand !=0
857 {
858 if(leftOperand_p->isFunction()) val=leftOperand_p->leftEval<T>(p, difop, d, m, nbw, np, extdata);
859 else val=leftOperand_p->leftEval<T>(difop,d,m,nbw);
860 if(rightOperand_p == 0) return;
861 if(rightOperand_p->isFunction()) val=rightOperand_p->rightEval<T>(p, val, d, m, nbw, np, extdata);
862 else val=rightOperand_p->rightEval<T>(val, d, m, nbw);
863 }
864 }
865
866 //evaluate operator from points p,q and shape values sv and its derivatives dsv
867 //(dimFun is the dimension of shape functions)
868 //case of a Kernel
869 template <typename T, typename K>
eval(const Point & p,const Point & q,const std::vector<K> & sv,const std::vector<std::vector<K>> & dsv,dimen_t dimFun,Vector<T> & val,dimen_t & d,dimen_t & m,const Vector<real_t> * nx,const Vector<real_t> * ny) const870 void OperatorOnUnknown::eval(const Point& p,const Point& q, const std::vector<K>& sv,
871 const std::vector< std::vector<K> >& dsv,
872 dimen_t dimFun, Vector<T>& val, dimen_t& d, dimen_t& m,
873 const Vector<real_t>* nx, const Vector<real_t>* ny) const
874 {
875 // compute differential operator
876 number_t nbw=sv.size()/dimFun;
877 d=dimFun;
878 Vector<K> difop;
879 difOp_p->eval(sv, dsv, d, m, difop, ny, coefs_);
880
881 //compute operator on unknown
882 if(leftOperand_p == 0 && rightOperand_p == 0) {val=difop; return;}
883 if(leftOperand_p != 0 && rightOperand_p == 0) {val=leftOperand_p->leftEval<T>(p,q,difop,d,m,nbw,nx,ny);return;}
884 if(leftOperand_p == 0 && rightOperand_p != 0) {val=rightOperand_p->rightEval<T>(p,q,difop,d,m,nbw,nx,ny);return;}
885 if(leftOperand_p->isFunction()) val=leftOperand_p->leftEval<T>(p,q,difop,d,m,nbw,nx,ny);
886 else val=leftOperand_p->leftEval<T>(difop,d,m,nbw);
887 if(rightOperand_p->isFunction()) val=rightOperand_p->rightEval<T>(p,q,val,d,m,nbw,nx,ny);
888 else val=rightOperand_p->rightEval<T>(val,d,m,nbw);
889 }
890
891 // print utility
892 std::ostream& operator<<(std::ostream&, const OperatorOnUnknown&); //!< outputs OperatorOnUnknown attributes
893
894 } // end of namespace xlifepp
895
896 #endif /* OPERATOR_ON_UNKNOWN_HPP */
897
898