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.cpp
19 \author E. Lunéville
20 \since 4 mar 2012
21 \date 9 may 2012
22
23 \brief Implementation of xlifepp::OperatorOnUnknown class members and related functions
24 */
25
26 //===============================================================================
27 // library dependencies
28 #include "OperatorOnUnknown.hpp"
29 //===============================================================================
30
31 namespace xlifepp
32 {
33
34 //------------------------------------------------------------------------------------
35 // OperatorOnUnknown member functions
36 //------------------------------------------------------------------------------------
37 // constructor of new OperatorOnUnknown
OperatorOnUnknown(const Unknown * un,DiffOpType ty)38 OperatorOnUnknown::OperatorOnUnknown(const Unknown* un, DiffOpType ty)
39 {
40 difOp_p = findDifferentialOperator(ty);
41 u_p = un;
42 conjugateUnknown_ = false;
43 if(un!=0 && un->conjugate())
44 {
45 conjugateUnknown_ = true;
46 un->conjugate(false); // reset the conjugate state flag of unknown
47 }
48 leftOperand_p = 0;
49 rightOperand_p = 0;
50 leftPriority_ = false;
51 dimsRes_ = dimPair(0,0);
52 setStructure();
53 }
54
OperatorOnUnknown(const Unknown & un,DiffOpType ty)55 OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, DiffOpType ty)
56 {
57 difOp_p = findDifferentialOperator(ty);
58 u_p = &un;
59 conjugateUnknown_ = false;
60 if(un.conjugate())
61 {
62 conjugateUnknown_ = true;
63 un.conjugate(false); // reset the conjugate state flag of unknown
64 }
65 leftOperand_p = 0;
66 rightOperand_p = 0;
67 leftPriority_ = false;
68 dimsRes_ = dimPair(0,0);
69 setStructure();
70 }
71
OperatorOnUnknown(const Unknown & un,DiffOpType ty,const std::vector<complex_t> & cs)72 OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, DiffOpType ty, const std::vector<complex_t>& cs)
73 {
74 difOp_p = findDifferentialOperator(ty);
75 coefs_ = cs;
76 u_p = &un;
77 conjugateUnknown_ = false;
78 if(un.conjugate())
79 {
80 conjugateUnknown_ = true;
81 un.conjugate(false); // reset the conjugate state flag of unknown
82 }
83 leftOperand_p = 0;
84 rightOperand_p = 0;
85 leftPriority_ = false;
86 dimsRes_ = dimPair(0,0);
87 setStructure();
88 }
89
90 // useful constructor : (fun op u) or (u op fun)
OperatorOnUnknown(const Unknown & un,const Function & fun,AlgebraicOperator aop,bool lef)91 OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, const Function& fun, AlgebraicOperator aop, bool lef)
92 {
93 difOp_p = findDifferentialOperator(_id);
94 u_p = &un;
95 conjugateUnknown_ = false;
96 if(un.conjugate())
97 {
98 conjugateUnknown_ = true;
99 un.conjugate(false); // reset the conjugate state flag of unknown
100 }
101 dimsRes_ = dimPair(0,0);
102 leftOperand_p = 0;
103 rightOperand_p = 0;
104 leftPriority_ = lef;
105 if(lef) { leftOperand_p = new Operand(fun, aop); }
106 else { rightOperand_p = new Operand(fun, aop); }
107 setStructure();
108 }
109
110 // useful constructor : (opfun op u) or (u op opfun)
OperatorOnUnknown(const Unknown & un,const OperatorOnFunction & opfun,AlgebraicOperator aop,bool lef)111 OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, const OperatorOnFunction& opfun, AlgebraicOperator aop, bool lef)
112 {
113 difOp_p = findDifferentialOperator(_id);
114 u_p = &un;
115 conjugateUnknown_ = false;
116 if(un.conjugate())
117 {
118 conjugateUnknown_ = true;
119 un.conjugate(false); // reset the conjugate state flag of unknown
120 }
121 dimsRes_ = dimPair(0,0);
122 leftOperand_p = 0;
123 rightOperand_p = 0;
124 leftPriority_ = lef;
125 if(lef) { leftOperand_p = new Operand(opfun, aop);}
126 else { rightOperand_p = new Operand(opfun, aop);}
127 setStructure();
128 }
129
130 // useful constructor : (val op u) or (u op val)
OperatorOnUnknown(const Unknown & un,const Value & val,AlgebraicOperator aop,bool lef)131 OperatorOnUnknown::OperatorOnUnknown(const Unknown& un, const Value& val, AlgebraicOperator aop, bool lef)
132 {
133 difOp_p = findDifferentialOperator(_id);
134 u_p = &un;
135 conjugateUnknown_ = false;
136 if(un.conjugate())
137 {
138 conjugateUnknown_ = true;
139 un.conjugate(false); // reset the conjugate state flag of unknown
140 }
141 dimsRes_ = dimPair(0,0);
142 leftOperand_p = 0;
143 rightOperand_p = 0;
144 leftPriority_ = lef;
145 if(lef) { leftOperand_p = new Operand(val, aop); }
146 else { rightOperand_p = new Operand(val, aop); }
147 setStructure();
148 }
149
150 // copy constructor
OperatorOnUnknown(const OperatorOnUnknown & opu)151 OperatorOnUnknown::OperatorOnUnknown(const OperatorOnUnknown& opu)
152 {
153 difOp_p = opu.difOp_p;
154 conjugateUnknown_ = opu.conjugateUnknown_;
155 u_p = opu.u_p;
156 type_ = opu.type_;
157 struct_ = opu.struct_;
158 dimsRes_ = opu.dimsRes_;
159 coefs_=opu.coefs_;
160 leftOperand_p = 0;
161 rightOperand_p = 0; // full copy
162 if(opu.leftOperand_p != 0) { leftOperand_p = new Operand(*opu.leftOperand_p); }
163 if(opu.rightOperand_p != 0) { rightOperand_p = new Operand(*opu.rightOperand_p); }
164 leftPriority_ = opu.leftPriority_;
165 }
166
167 // assignment operator
operator =(const OperatorOnUnknown & opu)168 OperatorOnUnknown& OperatorOnUnknown::operator=(const OperatorOnUnknown& opu)
169 {
170 difOp_p = opu.difOp_p;
171 conjugateUnknown_ = opu.conjugateUnknown_;
172 u_p = opu.u_p;
173 type_ = opu.type_;
174 struct_ = opu.struct_;
175 dimsRes_ = opu.dimsRes_;
176 leftPriority_ = opu.leftPriority_;
177 coefs_=opu.coefs_;
178 if(leftOperand_p != 0) { delete leftOperand_p; }
179 if(opu.leftOperand_p != 0) {leftOperand_p = new Operand(*opu.leftOperand_p);}
180 else {leftOperand_p =0;}
181 if(rightOperand_p != 0) { delete rightOperand_p; }
182 if(opu.rightOperand_p != 0) {rightOperand_p = new Operand(*opu.rightOperand_p);}
183 else {rightOperand_p =0;}
184 return *this;
185 }
186
187 // destructor
~OperatorOnUnknown()188 OperatorOnUnknown::~OperatorOnUnknown()
189 {
190 if(leftOperand_p != 0) { delete leftOperand_p; }
191 if(rightOperand_p != 0) { delete rightOperand_p; }
192 }
193
194 // set returned type, structure, dims of returned value and check consistency
195 // ### no longer used, replace by setStructure(), to be removed in next version ###
setReturnedType(const Unknown * un,DiffOpType ty)196 void OperatorOnUnknown::setReturnedType(const Unknown* un, DiffOpType ty)
197 {
198 type_ = _real;
199 struct_= _scalar;
200 dimsRes_ = dimPair(1,1);
201 if(un==0) return;
202
203 type_=un->space()->valueType(); //may be complex in case of spectral space with complex functions
204
205 dimen_t dimu = un->nbOfComponents(); // dimension of unknown
206 dimen_t dimf = un->space()->dimFun(); // dimension of shape functions
207 dimen_t dimv = std::max(dimu,dimf); // dimension for vector case
208 dimen_t dimd = un->space()->dimDomain(); // dimension of space domain
209 dimen_t dims = un->space()->dimPoint(); // dimension of points, may be greater than dimd
210
211 dimsRes_ = dimPair(dimu,1);
212
213 switch(ty)
214 {
215 case _grad:
216 case _gradS:
217 case _gradG:
218 {
219 struct_ = _vector;
220 dimsRes_=dimPair(dimd,1);
221 if(dimv > 1) { struct_ = _matrix; dimsRes_ = dimPair(dimd, dimv); }
222 break;
223 }
224 case _div:
225 case _divS:
226 case _divG:
227 {
228 struct_ = _scalar;
229 dimsRes_ =dimPair(1,1);
230 break;
231 }
232 case _curl:
233 case _curlS:
234 case _curlG:
235 {
236 struct_ = _vector;
237 dimsRes_ = dimPair(dimv,1);
238 break; // what about scalar curl in 2D
239 }
240 case _ntimes:
241 {
242 struct_ = _vector;
243 dimsRes_ = dimPair(dimv,1);
244 if(dimv > 1)
245 {
246 where("OperatorOnUnknown::setReturnedType(...)");
247 error("operator_baddiffop", words("diffop", ty));
248 }
249 break;
250 }
251 case _ndot:
252 {
253 struct_ = _scalar;
254 dimsRes_ = dimPair(1,1);
255 if(dimv != dims)
256 {
257 where("OperatorOnUnknown::setReturnedType(...)");
258 error("operator_baddiffop", words("diffop", ty));
259 }
260 break;
261 }
262 case _ncross:
263 case _ncrossncross:
264 case _ncrosscurl:
265 case _ndiv:
266 {
267 struct_ = _vector;
268 dimsRes_ = dimPair(dims,1);
269 break;
270 }
271 case _ndotgrad:
272 {
273 struct_ = _scalar;
274 dimsRes_ = dimPair(1,1);
275 if(dimv > 1) { struct_ = _vector; dimsRes_ = dimPair(dimv,1);}
276 break;
277 }
278 case _ncrossgrad:
279 {
280 struct_ = _scalar;
281 dimsRes_ = dimPair(1,1);
282 if(dims > 2) { struct_ = _vector; dimsRes_ = dimPair(dims,1);}
283 break;
284 }
285 case _epsilon:
286 case _epsilonG:
287 {
288 struct_ = _matrix;
289 dimsRes_ = dimPair(dimv,dimv);
290 break;
291 }
292 case _epsilonR:
293 {
294 struct_ = _vector;
295 dimsRes_ = dimPair((dimv*(dimv+1))/2,1);
296 break;
297 }
298 case _voigtToM:
299 {
300 struct_ = _matrix;
301 dimsRes_ = dimPair(dimv,dimv);
302 break;
303 }
304 default :
305 {
306 struct_ = _scalar;
307 dimsRes_ = dimPair(1,1);
308 if(dimv > 1) { struct_ = _vector; dimsRes_ = dimPair(dimv,1);}
309 }
310 }
311 }
312
313 // set returned type, structure, dims of returned value
setStructure()314 void OperatorOnUnknown::setStructure()
315 {
316 type_=_real;
317 struct_=_scalar;
318 dimsRes_=dimPair(0,0);
319 if(u_p==0) return; //no operator
320 // set value type
321 type_=u_p->space()->valueType(); //may be complex in case of spectral space with complex functions
322 if(type_==_real && leftOperand_p!=0) type_=leftOperand_p->valueType();
323 if(type_==_real && rightOperand_p!=0) type_=rightOperand_p->valueType();
324 //retry dimensions from unknown
325 dimen_t dimu = u_p->nbOfComponents(); // dimension of unknown
326 dimen_t dimf = u_p->space()->dimFun(); // dimension of shape functions
327 dimen_t dimv = std::max(dimu,dimf); // dimension for vector case
328 dimen_t dims = u_p->space()->dimPoint(); // dimension of points, may be greater than dimd
329 if(dimv > 1 && dimv < dims) dimv = dims; // shapevalues will be mapped from dimv to dims
330 Vector<real_t> np(1,1.);
331 if(dims==2) {np=Vector<real_t>(2,1.);}
332 else if(dims==3) {np=Vector<real_t>(3,1.);}
333 Point P=np, Q=P+1.;
334 dimen_t d=0, m=0;
335 //evaluate operator on a fake shape function
336 if(type_==_real)
337 {
338 std::vector<real_t> sv(dimv,1.);
339 Vector<real_t> val;
340 std::vector< std::vector<real_t> > dsv(std::max(dims,dimv),sv);
341 if(!hasFunction()) eval(sv, dsv, dimv, val, d, m,&np);
342 else if(hasFun())
343 {
344 if(normalRequired()) setNx(&np);
345 eval(P, sv, dsv, dimv, val, d, m, &np);
346 }
347 else if(hasKernel())
348 {
349 if(xnormalRequiredByOpKernel()) setNx(&np);
350 if(ynormalRequiredByOpKernel()) setNy(&np);
351 eval(P, Q, sv, dsv, dimv, val, d, m, &np, &np);
352 }
353 }
354 else // complex case
355 {
356 std::vector<real_t> sv(dimv,1.);
357 std::vector< std::vector<real_t> > dsv(std::max(dims,dimv),sv);
358 Vector<complex_t> val;
359 if(!hasFunction()) eval(sv, dsv, dimv, val, d, m,&np);
360 else if(hasFun())
361 {
362 if(normalRequired()) setNx(&np);
363 eval(P, sv, dsv, dimv, val, d, m, &np);
364 }
365 else if(hasKernel())
366 {
367 if(xnormalRequiredByOpKernel()) setNx(&np);
368 if(ynormalRequiredByOpKernel()) setNy(&np);
369 eval(P, Q, sv, dsv, dimv, val, d, m, &np, &np);
370 }
371 }
372 dimsRes_= dimPair(d/m,m);
373 if(dimsRes_.first==1) struct_=_scalar;
374 else if(dimsRes_.second==1) struct_=_vector;
375 else struct_=_matrix;
376 return;
377 }
378
379 //! order of polynom used by unknown (degree of basis)
unknownDegree() const380 number_t OperatorOnUnknown::unknownDegree() const
381 {
382 if(unknown()->space()->typeOfSpace()==_feSpace)
383 return unknown()->space()->feSpace()->maxDegree(); // new way, shape type no longer used
384 return 5; //not a FE space (no polynoms) (default order is chosen as 5)
385 }
386
387 /*! order of polynom used by unknown (degree od basis), decreased by order of derivation
388 as we use numtype (order of interpolation) we distinguish Pk or Qk using a ShapeType argument*/
degree() const389 number_t OperatorOnUnknown::degree() const
390 {
391 if(unknown()->space()->typeOfSpace()==_feSpace)
392 {
393 number_t d = unknown()->space()->feSpace()->maxDegree(); // degree of shape functions
394 if(hasFunction()) d*=2; //degree is double to take into account the function or kernel involved
395 d-=diffOrder(); //take into derivative
396 if(d==0) d=1; //min degree is 1
397 return d;
398 }
399 return 5; //not a FE space (no polynoms) (default order is chosen as 5)
400 }
401
402 //! order of differential operator involved
diffOrder() const403 dimen_t OperatorOnUnknown::diffOrder() const
404 {
405 if(difOp_p->type()==_epsilonG) return dimen_t(coefs_[0].real()); //special case
406 return difOp_p->order();
407 }
408
409 //! true if operator involves some functions or kernels
hasFunction() const410 bool OperatorOnUnknown::hasFunction() const
411 {
412 if(leftOperand_p != 0 && (leftOperand_p->isFunction() || leftOperand_p->isKernel())) { return true; }
413 if(rightOperand_p != 0 && (rightOperand_p->isFunction() || rightOperand_p->isKernel())) { return true; }
414 return false;
415 }
416
417 //! true if operator involves some functions or kernels
hasKernel() const418 bool OperatorOnUnknown::hasKernel() const
419 {
420 if(leftOperand_p != 0 && leftOperand_p->isKernel()) { return true; }
421 if(rightOperand_p != 0 && rightOperand_p->isKernel()) { return true; }
422 return false;
423 }
424
425 //! true if operator involves some functions
hasFun() const426 bool OperatorOnUnknown::hasFun() const
427 {
428 if(leftOperand_p != 0 && leftOperand_p->isFunction()) { return true; }
429 if(rightOperand_p != 0 && rightOperand_p->isFunction()) { return true; }
430 return false;
431 }
432
433 //! direct acces to Kernel (unique Function, left priority)
functionp() const434 const Function* OperatorOnUnknown::functionp() const
435 {
436 if(leftOperand_p != 0 && leftOperand_p->isFunction()) { return leftOperand_p->functionp(); }
437 if(rightOperand_p != 0 && rightOperand_p->isFunction()) { return rightOperand_p->functionp(); }
438 return 0;
439 }
440
441 //! direct acces to Kernel (unique Kernel, left priority)
kernelp() const442 const Kernel* OperatorOnUnknown::kernelp() const
443 {
444 if(leftOperand_p != 0 && leftOperand_p->isKernel()) { return leftOperand_p->kernelp(); }
445 if(rightOperand_p != 0 && rightOperand_p->isKernel()) { return rightOperand_p->kernelp(); }
446 return 0;
447 }
448
449 //! direct acces to Kernel (unique Kernel, left priority)
opkernelp() const450 const OperatorOnKernel* OperatorOnUnknown::opkernelp() const
451 {
452 if(leftOperand_p != 0 && leftOperand_p->isKernel()) { return leftOperand_p->opkernelp(); }
453 if(rightOperand_p != 0 && rightOperand_p->isKernel()) { return rightOperand_p->opkernelp(); }
454 return 0;
455 }
456
457
458 //! change Kernel in operator (unique Kernel, left priority)
changeKernel(Kernel * ker)459 void OperatorOnUnknown::changeKernel(Kernel* ker)
460 {
461 if(leftOperand_p != 0 && leftOperand_p->isKernel()) {leftOperand_p->changeKernel(ker); return;}
462 if(rightOperand_p != 0 && rightOperand_p->isKernel()) {rightOperand_p->changeKernel(ker); return;}
463 where("OperatorOnUnknown::changeKernel");
464 error("null_pointer","kernel");
465 }
466
467 //! true if operator involves left or right operand
hasOperand() const468 bool OperatorOnUnknown::hasOperand() const
469 {
470 return (leftOperand_p != 0 || rightOperand_p !=0);
471 }
472
473 //! true if normal is required in computation
normalRequired() const474 bool OperatorOnUnknown::normalRequired() const
475 {
476 if(difOp_p->normalRequired()) return true;
477 if(leftOperand_p != 0 && leftOperand_p->normalRequired()) return true;
478 if(rightOperand_p != 0 && rightOperand_p->normalRequired()) return true;
479 return false;
480 }
481
482 //! true if normal is required in computation
xnormalRequired() const483 bool OperatorOnUnknown::xnormalRequired() const
484 {
485 if(difOp_p->normalRequired()) return true;
486 if(leftOperand_p != 0 && leftOperand_p->xnormalRequired()) return true;
487 if(rightOperand_p != 0 && rightOperand_p->xnormalRequired()) return true;
488 return false;
489 }
490
491 //! true if normal is required in computation
ynormalRequired() const492 bool OperatorOnUnknown::ynormalRequired() const
493 {
494 if(difOp_p->normalRequired()) return true;
495 if(leftOperand_p != 0 && leftOperand_p->ynormalRequired()) return true;
496 if(rightOperand_p != 0 && rightOperand_p->ynormalRequired()) return true;
497 return false;
498 }
499
500 //! true if normal is required in computation by operator on Kernel
xnormalRequiredByOpKernel() const501 bool OperatorOnUnknown::xnormalRequiredByOpKernel() const
502 {
503 if(leftOperand_p != 0 && leftOperand_p->xnormalRequired()) return true;
504 if(rightOperand_p != 0 && rightOperand_p->xnormalRequired()) return true;
505 return false;
506 }
507
508 //! true if normal is required in computation by operator on Kernel
ynormalRequiredByOpKernel() const509 bool OperatorOnUnknown::ynormalRequiredByOpKernel() const
510 {
511 if(leftOperand_p != 0 && leftOperand_p->ynormalRequired()) return true;
512 if(rightOperand_p != 0 && rightOperand_p->ynormalRequired()) return true;
513 return false;
514 }
515
516 //!< true if GeomElemnt is required by operator
elementRequired() const517 bool OperatorOnUnknown::elementRequired() const
518 {
519 if(leftOperand_p != 0 && leftOperand_p->elementRequired()) return true;
520 if(rightOperand_p != 0 && rightOperand_p->elementRequired()) return true;
521 return false;
522 }
523
524 //! update value, structure and dimensions of result of (operator)(op)(valueType,strucType)
525 // dims : dimensions of value
updateReturnedType(AlgebraicOperator op,ValueType vt,StrucType st,dimPair dims,bool left)526 void OperatorOnUnknown::updateReturnedType(AlgebraicOperator op, ValueType vt, StrucType st, dimPair dims, bool left)
527 {
528 if(vt == _complex || type_ == _complex) { type_ = _complex; } //update value type
529
530 switch(op) // update structure type
531 {
532 case _product :
533 {
534 switch(struct_) // current opu struct
535 {
536 case _scalar :
537 {
538 struct_ = st; // scalar*any
539 return;
540 }
541 case _vector :
542 {
543 struct_ = _vector;
544 if(st == _scalar) // vector*scalar or scalar*vector
545 {
546 dimsRes_= dims;
547 return;
548 }
549 if(st == _matrix) // matrix*vector or vector*matrix
550 {
551 if(left) dimsRes_= dimPair(dims.first,1);
552 else dimsRes_= dimPair(dims.second,1);
553 return;
554 }
555 break;
556 }
557 case _matrix :
558 {
559 struct_ = _matrix;
560 if(st == _scalar) return; // scalar*matrix
561 if(st == _vector) // matrix*vector
562 {
563 if(left) dimsRes_= dimPair(dims.first,dimsRes_.second);
564 else dimsRes_= dimPair(dimsRes_.second, dims.second);
565 return;
566 }
567 if(st == _matrix) // matrix*matrix
568 {
569 if(left) dimsRes_= dimPair(dimsRes_.second,1);
570 else dimsRes_= dimPair(dimsRes_.first, 1);
571 return;
572 }
573 break;
574 }
575 default:
576 break;
577 }
578 break;
579 }
580 case _innerProduct :
581 {
582 switch(struct_) // current opu struct
583 {
584 case _scalar :
585 {
586 if(st == _scalar) // scalar.scalar
587 {
588 struct_ = _scalar;
589 return;
590 }
591 break;
592 }
593 case _vector :
594 {
595 if(st == _vector) // vector.vector
596 {
597 struct_ = _scalar;
598 dimsRes_=dimPair(1,1);
599 return;
600 }
601 break;
602 }
603 case _matrix :
604 break;
605 default:
606 break;
607 }
608 break;
609 }
610 case _crossProduct :
611 {
612 switch(struct_) // current opu struct
613 {
614 case _scalar :
615 case _matrix :
616 break;
617 case _vector :
618 {
619 if(st == _vector) // vector^vector
620 {
621 if(dims.first<=2) //return a scalar
622 {
623 struct_ = _scalar;
624 dimsRes_=dimPair(1,1);
625 return;
626 }
627 struct_ = _vector;
628 dimsRes_=dimPair(dims.first,1);
629 return;
630 }
631 }
632 default:
633 break;
634 }
635 }
636 break;
637 case _contractedProduct :
638 {
639 switch(struct_) // current opu struct
640 {
641 case _scalar :
642 case _vector :
643 break;
644 case _matrix :
645 {
646 if(st == _matrix) // matrix%matrix
647 //if matrix have the same size return a scalar (standard contracted product)
648 //else generalized contracted product returning a matrix of same size of original (see computation)
649 {
650 if(dimsRes_==dims)
651 {
652 struct_ = _scalar;
653 dimsRes_=dimPair(1,1);
654 return;
655 }
656 struct_ = _matrix;
657 return;
658 }
659 }
660 default:
661 break;
662 }
663 }
664 }
665 if(left)
666 {
667 where("OperatorOnUnknown::updateReturnedType(...)");
668 error("operator_badop", words("structure", struct_), words("algop", op), words("structure", st));
669 }
670 else
671 {
672 where("OperatorOnUnknown::updateReturnedType(...)");
673 error("operator_badop", words("structure", st), words("algop", op), words("structure", struct_));
674 }
675 }
676
677 //------------------------------------------------------------------------------------
678 // external functions relates to OperatorOnUnknown
679 //------------------------------------------------------------------------------------
680
681 // user interface functions
682 // be cautious with memory, OperatorOnUnknown are created on dynamic memory
683 // at the end of the chaining construction, be sure to free it safely
684
685 //==========================================================================================
686 // operators involving DifferentialOperator and Unknown (create new OperatorOnUnknown)
687 //==========================================================================================
id(const Unknown & un)688 OperatorOnUnknown& id(const Unknown& un)
689 {
690 return*(new OperatorOnUnknown(&un, _id));
691 }
d0(const Unknown & un)692 OperatorOnUnknown& d0(const Unknown& un)
693 {
694 return*(new OperatorOnUnknown(&un, _d0));
695 }
dt(const Unknown & un)696 OperatorOnUnknown& dt(const Unknown& un)
697 {
698 return*(new OperatorOnUnknown(&un, _d0));
699 }
d1(const Unknown & un)700 OperatorOnUnknown& d1(const Unknown& un)
701 {
702 return*(new OperatorOnUnknown(&un, _d1));
703 }
dx(const Unknown & un)704 OperatorOnUnknown& dx(const Unknown& un)
705 {
706 return*(new OperatorOnUnknown(&un, _d1));
707 }
d2(const Unknown & un)708 OperatorOnUnknown& d2(const Unknown& un)
709 {
710 return*(new OperatorOnUnknown(&un, _d2));
711 }
dy(const Unknown & un)712 OperatorOnUnknown& dy(const Unknown& un)
713 {
714 return*(new OperatorOnUnknown(&un, _d2));
715 }
d3(const Unknown & un)716 OperatorOnUnknown& d3(const Unknown& un)
717 {
718 return*(new OperatorOnUnknown(&un, _d3));
719 }
dz(const Unknown & un)720 OperatorOnUnknown& dz(const Unknown& un)
721 {
722 return*(new OperatorOnUnknown(&un, _d3));
723 }
grad(const Unknown & un)724 OperatorOnUnknown& grad(const Unknown& un)
725 {
726 return *(new OperatorOnUnknown(&un, _grad));
727 }
nabla(const Unknown & un)728 OperatorOnUnknown& nabla(const Unknown& un)
729 {
730 return*(new OperatorOnUnknown(&un, _grad));
731 }
div(const Unknown & un)732 OperatorOnUnknown& div(const Unknown& un)
733 {
734 return*(new OperatorOnUnknown(&un, _div));
735 }
curl(const Unknown & un)736 OperatorOnUnknown& curl(const Unknown& un)
737 {
738 return*(new OperatorOnUnknown(&un, _curl));
739 }
rot(const Unknown & un)740 OperatorOnUnknown& rot(const Unknown& un)
741 {
742 return*(new OperatorOnUnknown(&un, _curl));
743 }
gradS(const Unknown & un)744 OperatorOnUnknown& gradS(const Unknown& un)
745 {
746 return*(new OperatorOnUnknown(&un, _gradS));
747 }
nablaS(const Unknown & un)748 OperatorOnUnknown& nablaS(const Unknown& un)
749 {
750 return*(new OperatorOnUnknown(&un, _gradS));
751 }
divS(const Unknown & un)752 OperatorOnUnknown& divS(const Unknown& un)
753 {
754 return*(new OperatorOnUnknown(&un, _divS));
755 }
curlS(const Unknown & un)756 OperatorOnUnknown& curlS(const Unknown& un)
757 {
758 return*(new OperatorOnUnknown(&un, _curlS));
759 }
rotS(const Unknown & un)760 OperatorOnUnknown& rotS(const Unknown& un)
761 {
762 return*(new OperatorOnUnknown(&un, _curlS));
763 }
nx(const Unknown & un)764 OperatorOnUnknown& nx(const Unknown& un)
765 {
766 return*(new OperatorOnUnknown(&un, _ntimes));
767 }
ndot(const Unknown & un)768 OperatorOnUnknown& ndot(const Unknown& un)
769 {
770 return*(new OperatorOnUnknown(&un, _ndot));
771 }
ndotgrad(const Unknown & un)772 OperatorOnUnknown& ndotgrad(const Unknown& un)
773 {
774 return*(new OperatorOnUnknown(&un, _ndotgrad));
775 }
ncrossgrad(const Unknown & un)776 OperatorOnUnknown& ncrossgrad(const Unknown& un)
777 {
778 return*(new OperatorOnUnknown(&un, _ncrossgrad));
779 }
ncross(const Unknown & un)780 OperatorOnUnknown& ncross(const Unknown& un)
781 {
782 return*(new OperatorOnUnknown(&un, _ncross));
783 }
ncrossncross(const Unknown & un)784 OperatorOnUnknown& ncrossncross(const Unknown& un)
785 {
786 return*(new OperatorOnUnknown(&un, _ncrossncross));
787 }
ncrossntimes(const Unknown & un)788 OperatorOnUnknown& ncrossntimes(const Unknown& un)
789 {
790 return*(new OperatorOnUnknown(&un, _ncrossntimes));
791 }
792 //OperatorOnUnknown& ncrossndot(const Unknown& un)
793 //{
794 // return*(new OperatorOnUnknown(&un, _ncrossndot));
795 //}
ntimesndot(const Unknown & un)796 OperatorOnUnknown& ntimesndot(const Unknown& un)
797 {
798 return*(new OperatorOnUnknown(&un, _ntimesndot));
799 }
ndiv(const Unknown & un)800 OperatorOnUnknown& ndiv(const Unknown& un)
801 {
802 return*(new OperatorOnUnknown(&un, _ndiv));
803 }
ncrosscurl(const Unknown & un)804 OperatorOnUnknown& ncrosscurl(const Unknown& un)
805 {
806 return*(new OperatorOnUnknown(&un, _ncrosscurl));
807 }
epsilon(const Unknown & un)808 OperatorOnUnknown& epsilon(const Unknown& un)
809 {
810 return*(new OperatorOnUnknown(&un, _epsilon));
811 }
812
epsilonR(const Unknown & un)813 OperatorOnUnknown& epsilonR(const Unknown& un)
814 {
815 return*(new OperatorOnUnknown(&un, _epsilonR));
816 }
817
voigtToM(const Unknown & un)818 OperatorOnUnknown& voigtToM(const Unknown& un) //special operator : convert Voigt vector to tensor matrix
819 {
820 return*(new OperatorOnUnknown(&un, _voigtToM));
821 }
822
823 // generalized differential operator
gradG(const Unknown & un,const complex_t & ax,const complex_t & ay,const complex_t & az,const complex_t & at)824 OperatorOnUnknown& gradG(const Unknown& un, const complex_t& ax, const complex_t& ay,
825 const complex_t& az, const complex_t& at)
826 {
827 std::vector<complex_t> cs(4);
828 cs[0]=ax; cs[1]=ay; cs[2]=az; cs[3]=at;
829 OperatorOnUnknown* opu = new OperatorOnUnknown(un, _gradG, cs);
830 return *opu;
831 }
832
nablaG(const Unknown & un,const complex_t & ax,const complex_t & ay,const complex_t & az,const complex_t & at)833 OperatorOnUnknown& nablaG(const Unknown& un, const complex_t& ax, const complex_t& ay,
834 const complex_t& az, const complex_t& at)
835 {
836 return gradG(un, ax, ay, az, at);
837 }
838
divG(const Unknown & un,const complex_t & ax,const complex_t & ay,const complex_t & az,const complex_t & at)839 OperatorOnUnknown& divG(const Unknown& un, const complex_t& ax, const complex_t& ay,
840 const complex_t& az, const complex_t& at)
841 {
842 std::vector<complex_t> cs(4);
843 cs[0]=ax; cs[1]=ay; cs[2]=az; cs[3]=at;
844 OperatorOnUnknown* opu = new OperatorOnUnknown(un, _divG, cs);
845 return *opu;
846 }
847
curlG(const Unknown & un,const complex_t & ax,const complex_t & ay,const complex_t & az,const complex_t & at)848 OperatorOnUnknown& curlG(const Unknown& un, const complex_t& ax, const complex_t& ay,
849 const complex_t& az, const complex_t& at)
850 {
851 std::vector<complex_t> cs(4);
852 cs[0]=ax; cs[1]=ay; cs[2]=az; cs[3]=at;
853 OperatorOnUnknown* opu = new OperatorOnUnknown(un, _curlG, cs);
854 return *opu;
855 }
856
rotG(const Unknown & un,const complex_t & ax,const complex_t & ay,const complex_t & az,const complex_t & at)857 OperatorOnUnknown& rotG(const Unknown& un, const complex_t& ax, const complex_t& ay,
858 const complex_t& az, const complex_t& at)
859 {
860 return curlG(un, ax, ay, az, at);
861 }
862
epsilonG(const Unknown & un,const complex_t & id,const complex_t & ax,const complex_t & ay,const complex_t & az)863 OperatorOnUnknown& epsilonG(const Unknown& un, const complex_t& id, const complex_t& ax,
864 const complex_t& ay, const complex_t& az)
865 {
866 std::vector<complex_t> cs(4);
867 cs[0]=id; cs[1]=ax; cs[2]=ay; cs[3]=az;
868 OperatorOnUnknown* opu = new OperatorOnUnknown(un, _epsilonG, cs);
869 return *opu;
870 }
871
mean(const Unknown & un)872 OperatorOnUnknown& mean(const Unknown& un)
873 {
874 return*(new OperatorOnUnknown(&un, _mean));
875 }
876
jump(const Unknown & un)877 OperatorOnUnknown& jump(const Unknown& un)
878 {
879 return*(new OperatorOnUnknown(&un, _jump));
880 }
881
882 //==========================================================================================
883 // operators involving UnitaryVector and Unknown (create new OperatorOnUnknown)
884 //==========================================================================================
operator *(UnitaryVector v,const Unknown & un)885 OperatorOnUnknown& operator*(UnitaryVector v, const Unknown& un) // n*u same as nx(u)
886 {
887 if(v == _n) return *(new OperatorOnUnknown(&un, _ntimes));
888 if(v == _ncrossn) return *(new OperatorOnUnknown(&un, _ncrossntimes));
889 where("UnitaryVector * Unknown");
890 error("operator_unexpected", words("diffop",v)+ "? * u");
891 return *(new OperatorOnUnknown()); // dummy return
892 }
893
operator |(UnitaryVector v,const Unknown & un)894 OperatorOnUnknown& operator|(UnitaryVector v, const Unknown& un) // n|u same as ndot(u)
895 {
896 if(v == _n) return *(new OperatorOnUnknown(&un, _ndot));
897 // if(v == _ncrossn) return *(new OperatorOnUnknown(&un, _ncrossndot));
898 where("UnitaryVector | Unknown");
899 error("operator_unexpected", words("diffop",v)+ "? | u");
900 return *(new OperatorOnUnknown()); // dummy return
901 }
902
operator ^(UnitaryVector v,const Unknown & un)903 OperatorOnUnknown& operator^(UnitaryVector v, const Unknown& un) // n^u same as ncross(u)
904 {
905 if(v == _n) return *(new OperatorOnUnknown(&un, _ncross));
906 if(v == _ncrossn) return *(new OperatorOnUnknown(&un, _ncrossncross));
907 where("UnitaryVector ^ Unknown");
908 error("operator_unexpected", words("diffop",v)+ "? ^ u");
909 return *(new OperatorOnUnknown()); // dummy return
910 }
911
operator ^(const Unknown & un,UnitaryVector v)912 OperatorOnUnknown& operator^(const Unknown& un, UnitaryVector v) // u^n same as -ncross(u)
913 {
914 OperatorOnUnknown* opu=0;
915 if(v == _n) opu = new OperatorOnUnknown(&un, _ncross);
916 if(v == _ncrossn) opu = new OperatorOnUnknown(&un, _ncrossncross);
917 where("Unknown ^ UnitaryVector");
918 if(opu==0) error("operator_unexpected", words("diffop",v)+ "u ^ ?");
919 opu->leftOperand() = new Operand(-1.,_product);
920 return *opu;
921 }
922
operator *(const Unknown & un,UnitaryVector v)923 OperatorOnUnknown& operator*(const Unknown& un, UnitaryVector v) // u*n same as nx(u)
924 {return v * un;}
925
operator |(const Unknown & un,UnitaryVector v)926 OperatorOnUnknown& operator|(const Unknown& un, UnitaryVector v) // u|n same as ndot(u)
927 {return v | un;}
928
operator |(UnitaryVector v,OperatorOnUnknown & opu)929 OperatorOnUnknown& operator|(UnitaryVector v, OperatorOnUnknown& opu) // n|grad(u) same as ndotgrad(u)
930 {
931 if(v == _n && opu.difOpType() == _grad) {opu.difOp_() = findDifferentialOperator(_ndotgrad); return opu;}
932 if(v == _n && opu.difOpType() == _id) {opu.difOp_() = findDifferentialOperator(_ndot); return opu;}
933 // if(v == _ncrossn && opu.difOpType() == _id) {opu.difOp_() = findDifferentialOperator(_ncrossndot); return opu;}
934 where("UnitaryVector | OperatorOnUnknown");
935 error("operator_unexpected", words("diffop",v)+ " | "+ words("diffop",opu.difOpType()));
936 return opu;
937 }
938
operator *(UnitaryVector v,OperatorOnUnknown & opu)939 OperatorOnUnknown& operator*(UnitaryVector v, OperatorOnUnknown& opu) // n*, n*div, n^n*, n*n|
940 {
941 if(v == _n && opu.difOpType() == _id) {opu.difOp_() = findDifferentialOperator(_ntimes); return opu;}
942 if(v == _n && opu.difOpType() == _div) {opu.difOp_() = findDifferentialOperator(_ndiv); return opu;}
943 if(v == _ncrossn && opu.difOpType() == _id) {opu.difOp_() = findDifferentialOperator(_ncrossntimes); return opu;}
944 if(v == _n && opu.difOpType()==_ndot) {opu.difOp_() = findDifferentialOperator(_ntimesndot); return opu;}
945 where("UnitaryVector * OperatorOnUnknown");
946 error("operator_unexpected", words("diffop",v)+ " * "+ words("diffop",opu.difOpType()));
947 return opu;
948 }
949
operator ^(UnitaryVector v,OperatorOnUnknown & opu)950 OperatorOnUnknown& operator^(UnitaryVector v, OperatorOnUnknown& opu) // n^(n^u) or n^curl(u) or n^grad(u) same as ncrossncross(u) or ncrosscurl(u) or ncrossgrad(u)
951 {
952 DiffOpType ty = opu.difOpType();
953 if(v != _n)
954 {
955 where("UnitaryVector ^ OperatorOnUnknown");
956 error("operator_unexpected", words("normal",v)+ "? ^ "+ words("diffop",ty));
957 }
958 switch(ty)
959 {
960 case _ncross: opu.difOp_() = findDifferentialOperator(_ncrossncross); break;
961 case _grad: opu.difOp_() = findDifferentialOperator(_ncrossgrad); break;
962 case _curl: opu.difOp_() = findDifferentialOperator(_ncrosscurl); break;
963 default:
964 where("UnitaryVector ^ OperatorOnUnknown");
965 error("operator_unexpected", words("normal",v)+ " ^ "+ words("diffop",ty)+"?");
966 }
967 return opu;
968 }
969
operator ^(OperatorOnUnknown & opu,UnitaryVector v)970 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, UnitaryVector v) // (n^u)^n or curl(u)^n or grad(u)^n same as -ncrossncross(u) or -ncrosscurl(u) or -ncrossgrad(u)
971 {
972 DiffOpType ty = opu.difOpType();
973 if(v != _n)
974 {
975 where("OperatorOnUnknown ^ UnitaryVector");
976 error("operator_unexpected", words("diffop",ty)+ " ^ " + words("normal",v) + "?");
977 }
978 switch(ty)
979 {
980 case _ncross: opu.difOp_() = findDifferentialOperator(_ncrossncross); break;
981 case _grad: opu.difOp_() = findDifferentialOperator(_ncrossgrad); break;
982 case _curl: opu.difOp_() = findDifferentialOperator(_ncrosscurl); break;
983 default:
984 where("OperatorOnUnknown ^ UnitaryVector");
985 error("operator_unexpected", words("diffop",ty)+ " ^ "+words("normal",v)+ "?");
986 }
987 if(opu.leftOperand()==0) opu.leftOperand()= new Operand(-1.,_product);
988 else
989 {
990 if(opu.leftOperand()->isValue()) opu.leftOperand()->value().opposite();
991 else
992 {
993 where("OperatorOnUnknown ^ UnitaryVector");
994 error("opposite_operand_impossible");
995 }
996 }
997 return opu;
998 }
999
operator |(OperatorOnUnknown & opu,UnitaryVector v)1000 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, UnitaryVector v) // only grad(u)|n same as ndotgrad(u)
1001 {return v | opu;}
1002
operator *(OperatorOnUnknown & opu,UnitaryVector v)1003 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, UnitaryVector v) // only div*n same as ndiv(u)
1004 {return v * opu;}
1005
1006 //==========================================================================================
1007 // operators involving Unknown and Function (create new OperatorOnUnknown)
1008 //==========================================================================================
operator *(const Unknown & un,const Function & fun)1009 OperatorOnUnknown& operator*(const Unknown& un, const Function& fun) // u*F
1010 {return *(new OperatorOnUnknown(un, fun, _product, false));}
1011
operator |(const Unknown & un,const Function & fun)1012 OperatorOnUnknown& operator|(const Unknown& un, const Function& fun) // u|F
1013 {return *(new OperatorOnUnknown(un, fun, _innerProduct, false));}
1014
operator ^(const Unknown & un,const Function & fun)1015 OperatorOnUnknown& operator^(const Unknown& un, const Function& fun) // u^F
1016 {return *(new OperatorOnUnknown(un, fun, _crossProduct, false));}
1017
operator %(const Unknown & un,const Function & fun)1018 OperatorOnUnknown& operator%(const Unknown& un, const Function& fun) // u%F
1019 {return *(new OperatorOnUnknown(un, fun, _contractedProduct, false));}
1020
operator *(const Function & fun,const Unknown & un)1021 OperatorOnUnknown& operator*(const Function& fun, const Unknown& un) // F*u
1022 {return *(new OperatorOnUnknown(un, fun, _product, true));}
1023
operator |(const Function & fun,const Unknown & un)1024 OperatorOnUnknown& operator|(const Function& fun, const Unknown& un) // F|u
1025 {return *(new OperatorOnUnknown(un, fun, _innerProduct, true));}
1026
operator ^(const Function & fun,const Unknown & un)1027 OperatorOnUnknown& operator^(const Function& fun, const Unknown& un) // F^u
1028 {return *(new OperatorOnUnknown(un, fun, _crossProduct, true));}
1029
operator %(const Function & fun,const Unknown & un)1030 OperatorOnUnknown& operator%(const Function& fun, const Unknown& un) // F%u
1031 {return *(new OperatorOnUnknown(un, fun, _contractedProduct, true));}
1032
1033 //==========================================================================================
1034 // operators involving Unknown and OperatorOnFunction (create new OperatorOnUnknown)
1035 //==========================================================================================
operator *(const Unknown & un,const OperatorOnFunction & opfun)1036 OperatorOnUnknown& operator*(const Unknown& un, const OperatorOnFunction& opfun) // u*op(F)
1037 {return *(new OperatorOnUnknown(un, opfun, _product, false));}
1038
operator |(const Unknown & un,const OperatorOnFunction & opfun)1039 OperatorOnUnknown& operator|(const Unknown& un, const OperatorOnFunction& opfun) // u|op(F)
1040 {return *(new OperatorOnUnknown(un, opfun, _innerProduct, false));}
1041
operator ^(const Unknown & un,const OperatorOnFunction & opfun)1042 OperatorOnUnknown& operator^(const Unknown& un, const OperatorOnFunction& opfun) // u^op(F)
1043 {return *(new OperatorOnUnknown(un, opfun, _crossProduct, false));}
1044
operator %(const Unknown & un,const OperatorOnFunction & opfun)1045 OperatorOnUnknown& operator%(const Unknown& un, const OperatorOnFunction& opfun) // u%op(F)
1046 {return *(new OperatorOnUnknown(un, opfun, _contractedProduct, false));}
1047
operator *(const OperatorOnFunction & opfun,const Unknown & un)1048 OperatorOnUnknown& operator*(const OperatorOnFunction& opfun, const Unknown& un) // op(F)*u
1049 {return *(new OperatorOnUnknown(un, opfun, _product, true));}
1050
operator |(const OperatorOnFunction & opfun,const Unknown & un)1051 OperatorOnUnknown& operator|(const OperatorOnFunction& opfun, const Unknown& un) // op(F)|u
1052 {return *(new OperatorOnUnknown(un, opfun, _innerProduct, true));}
1053
operator ^(const OperatorOnFunction & opfun,const Unknown & un)1054 OperatorOnUnknown& operator^(const OperatorOnFunction& opfun, const Unknown& un) // op(F)^u
1055 {return *(new OperatorOnUnknown(un, opfun, _crossProduct, true));}
1056
operator %(const OperatorOnFunction & opfun,const Unknown & un)1057 OperatorOnUnknown& operator%(const OperatorOnFunction& opfun, const Unknown& un) // op(F)%u
1058 {return *(new OperatorOnUnknown(un, opfun, _contractedProduct, true));}
1059
1060 //==========================================================================================
1061 // operators involving Unknown and Value (create new OperatorOnUnknown)
1062 //==========================================================================================
operator *(const Unknown & un,const Value & val)1063 OperatorOnUnknown& operator*(const Unknown& un, const Value& val) // u*val
1064 {return *(new OperatorOnUnknown(un, val, _product, false));}
1065
operator |(const Unknown & un,const Value & val)1066 OperatorOnUnknown& operator|(const Unknown& un, const Value& val) // u|val
1067 {return *(new OperatorOnUnknown(un, val, _innerProduct, false));}
1068
operator ^(const Unknown & un,const Value & val)1069 OperatorOnUnknown& operator^(const Unknown& un, const Value& val) // u^val
1070 {return *(new OperatorOnUnknown(un, val, _crossProduct, false));}
1071
operator %(const Unknown & un,const Value & val)1072 OperatorOnUnknown& operator%(const Unknown& un, const Value& val) // u%val
1073 {return *(new OperatorOnUnknown(un, val, _contractedProduct, false));}
1074
operator *(const Value & val,const Unknown & un)1075 OperatorOnUnknown& operator*(const Value& val, const Unknown& un) // val*u
1076 {return *(new OperatorOnUnknown(un, val, _product, true));}
1077
operator |(const Value & val,const Unknown & un)1078 OperatorOnUnknown& operator|(const Value& val, const Unknown& un) // val|u
1079 {return *(new OperatorOnUnknown(un, val, _innerProduct, true));}
1080
operator ^(const Value & val,const Unknown & un)1081 OperatorOnUnknown& operator^(const Value& val, const Unknown& un) // val^u
1082 {return *(new OperatorOnUnknown(un, val, _crossProduct, true));}
1083
operator %(const Value & val,const Unknown & un)1084 OperatorOnUnknown& operator%(const Value& val, const Unknown& un) // val%u
1085 {return *(new OperatorOnUnknown(un, val, _contractedProduct, true));}
1086
1087 //==========================================================================================
1088 // operators involving Unknown and explicit value (create new OperatorOnUnknown)
1089 //==========================================================================================
operator *(const Unknown & un,const real_t & val)1090 OperatorOnUnknown& operator*(const Unknown& un, const real_t& val)
1091 {return fromUnknownVal(un,val,_product);}
1092
operator |(const Unknown & un,const real_t & val)1093 OperatorOnUnknown& operator|(const Unknown& un, const real_t& val)
1094 {return fromUnknownVal(un,val,_innerProduct);}
1095
operator ^(const Unknown & un,const real_t & val)1096 OperatorOnUnknown& operator^(const Unknown& un, const real_t& val)
1097 {return fromUnknownVal(un,val,_crossProduct);}
1098
operator %(const Unknown & un,const real_t & val)1099 OperatorOnUnknown& operator%(const Unknown& un, const real_t& val)
1100 {return fromUnknownVal(un,val,_contractedProduct);}
1101
operator *(const real_t & val,const Unknown & un)1102 OperatorOnUnknown& operator*(const real_t& val, const Unknown& un)
1103 {return fromValUnknown(un,val,_product);}
1104
operator |(const real_t & val,const Unknown & un)1105 OperatorOnUnknown& operator|(const real_t& val, const Unknown& un)
1106 {return fromValUnknown(un,val,_innerProduct);}
1107
operator ^(const real_t & val,const Unknown & un)1108 OperatorOnUnknown& operator^(const real_t& val, const Unknown& un)
1109 {return fromValUnknown(un,val,_crossProduct);}
1110
operator %(const real_t & val,const Unknown & un)1111 OperatorOnUnknown& operator%(const real_t& val, const Unknown& un)
1112 {return fromValUnknown(un,val,_contractedProduct);}
1113
operator *(const Unknown & un,const complex_t & val)1114 OperatorOnUnknown& operator*(const Unknown& un, const complex_t& val)
1115 {return fromUnknownVal(un,val,_product);}
1116
operator |(const Unknown & un,const complex_t & val)1117 OperatorOnUnknown& operator|(const Unknown& un, const complex_t& val)
1118 {return fromUnknownVal(un,val,_innerProduct);}
1119
operator ^(const Unknown & un,const complex_t & val)1120 OperatorOnUnknown& operator^(const Unknown& un, const complex_t& val)
1121 {return fromUnknownVal(un,val,_crossProduct);}
1122
operator %(const Unknown & un,const complex_t & val)1123 OperatorOnUnknown& operator%(const Unknown& un, const complex_t& val)
1124 {return fromUnknownVal(un,val,_contractedProduct);}
1125
operator *(const complex_t & val,const Unknown & un)1126 OperatorOnUnknown& operator*(const complex_t& val, const Unknown& un)
1127 {return fromValUnknown(un,val,_product);}
1128
operator |(const complex_t & val,const Unknown & un)1129 OperatorOnUnknown& operator|(const complex_t& val, const Unknown& un)
1130 {return fromValUnknown(un,val,_innerProduct);}
1131
operator ^(const complex_t & val,const Unknown & un)1132 OperatorOnUnknown& operator^(const complex_t& val, const Unknown& un)
1133 {return fromValUnknown(un,val,_crossProduct);}
1134
operator %(const complex_t & val,const Unknown & un)1135 OperatorOnUnknown& operator%(const complex_t& val, const Unknown& un)
1136 {return fromValUnknown(un,val,_contractedProduct);}
1137
1138
1139 //==========================================================================================
1140 // operators involving OperatorOnUnknown and Function (update current OperatorOnUnknown)
1141 //==========================================================================================
1142 // update left and right operation with Function : op(u) aop Function or Function aop op(u)
updateRight(OperatorOnUnknown & opu,const Function & fun,AlgebraicOperator aop)1143 OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const Function& fun, AlgebraicOperator aop)
1144 {
1145 if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
1146 //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), false);
1147 opu.rightOperand() = new Operand(fun, aop);
1148 if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
1149 opu.setStructure();
1150 return opu;
1151 }
updateLeft(OperatorOnUnknown & opu,const Function & fun,AlgebraicOperator aop)1152 OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const Function& fun, AlgebraicOperator aop)
1153 {
1154 if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
1155 //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), true);
1156 opu.leftOperand() = new Operand(fun, aop);
1157 if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
1158 opu.setStructure();
1159 return opu;
1160 }
1161
1162 // left and right operation with Function : op(u) aop F or F aop op(u)
operator *(OperatorOnUnknown & opu,const Function & fun)1163 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const Function& fun) // d(u) * F
1164 {
1165 return updateRight(opu, fun, _product);
1166 }
operator |(OperatorOnUnknown & opu,const Function & fun)1167 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const Function& fun) // d(u) | F
1168 {
1169 return updateRight(opu, fun, _innerProduct);
1170 }
operator ^(OperatorOnUnknown & opu,const Function & fun)1171 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const Function& fun) // d(u) ^ F
1172 {
1173 return updateRight(opu, fun, _crossProduct);
1174 }
operator %(OperatorOnUnknown & opu,const Function & fun)1175 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const Function& fun) // d(u) % F
1176 {
1177 return updateRight(opu, fun, _contractedProduct);
1178 }
operator *(const Function & fun,OperatorOnUnknown & opu)1179 OperatorOnUnknown& operator*(const Function& fun, OperatorOnUnknown& opu) // F * d(u)
1180 {
1181 return updateLeft(opu, fun, _product);
1182 }
operator |(const Function & fun,OperatorOnUnknown & opu)1183 OperatorOnUnknown& operator|(const Function& fun, OperatorOnUnknown& opu) // F | d(u)
1184 {
1185 return updateLeft(opu, fun, _innerProduct);
1186 }
operator ^(const Function & fun,OperatorOnUnknown & opu)1187 OperatorOnUnknown& operator^(const Function& fun, OperatorOnUnknown& opu) // F ^ d(u)
1188 {
1189 return updateLeft(opu, fun, _crossProduct);
1190 }
operator %(const Function & fun,OperatorOnUnknown & opu)1191 OperatorOnUnknown& operator%(const Function& fun, OperatorOnUnknown& opu) // F % d(u)
1192 {
1193 return updateLeft(opu, fun, _contractedProduct);
1194 }
1195
1196 //==========================================================================================
1197 // operators involving OperatorOnUnknown and OperatorOnFunction (update current OperatorOnUnknown)
1198 //==========================================================================================
1199 // update left and right operation with Function : op(u) aop op(F) or op(F) aop op(u)
updateRight(OperatorOnUnknown & opu,const OperatorOnFunction & fun,AlgebraicOperator aop)1200 OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const OperatorOnFunction& fun, AlgebraicOperator aop)
1201 {
1202 if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
1203 //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), false);
1204 opu.rightOperand() = new Operand(fun, aop);
1205 if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
1206 opu.setStructure();
1207 return opu;
1208 }
updateLeft(OperatorOnUnknown & opu,const OperatorOnFunction & fun,AlgebraicOperator aop)1209 OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const OperatorOnFunction& fun, AlgebraicOperator aop)
1210 {
1211 if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
1212 //opu.updateReturnedType(aop, fun.valueType(), fun.strucType(), fun.dims(), true);
1213 opu.leftOperand() = new Operand(fun, aop);
1214 if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
1215 opu.setStructure();
1216 return opu;
1217 }
1218
1219 // left and right operation with OperatorOnFunction : op(u) aop F or F aop op(u)
operator *(OperatorOnUnknown & opu,const OperatorOnFunction & fun)1220 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) * F
1221 {
1222 return updateRight(opu, fun, _product);
1223 }
operator |(OperatorOnUnknown & opu,const OperatorOnFunction & fun)1224 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) | F
1225 {
1226 return updateRight(opu, fun, _innerProduct);
1227 }
operator ^(OperatorOnUnknown & opu,const OperatorOnFunction & fun)1228 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) ^ F
1229 {
1230 return updateRight(opu, fun, _crossProduct);
1231 }
operator %(OperatorOnUnknown & opu,const OperatorOnFunction & fun)1232 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const OperatorOnFunction& fun) // d(u) % F
1233 {
1234 return updateRight(opu, fun, _contractedProduct);
1235 }
operator *(const OperatorOnFunction & fun,OperatorOnUnknown & opu)1236 OperatorOnUnknown& operator*(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F * d(u)
1237 {
1238 return updateLeft(opu, fun, _product);
1239 }
operator |(const OperatorOnFunction & fun,OperatorOnUnknown & opu)1240 OperatorOnUnknown& operator|(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F | d(u)
1241 {
1242 return updateLeft(opu, fun, _innerProduct);
1243 }
operator ^(const OperatorOnFunction & fun,OperatorOnUnknown & opu)1244 OperatorOnUnknown& operator^(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F ^ d(u)
1245 {
1246 return updateLeft(opu, fun, _crossProduct);
1247 }
operator %(const OperatorOnFunction & fun,OperatorOnUnknown & opu)1248 OperatorOnUnknown& operator%(const OperatorOnFunction& fun, OperatorOnUnknown& opu) // F % d(u)
1249 {
1250 return updateLeft(opu, fun, _contractedProduct);
1251 }
1252
1253 //==========================================================================================
1254 // operators involving OperatorOnUnknown and scalar (update current OperatorOnUnknown)
1255 //==========================================================================================
1256 // left and right operation with real scalar : op(u) aop r or r aop op(u)
operator *(OperatorOnUnknown & opu,const real_t & val)1257 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const real_t& val)
1258 {return updateRight(opu, val, _product);}
1259
operator |(OperatorOnUnknown & opu,const real_t & val)1260 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const real_t& val)
1261 {return updateRight(opu, val, _innerProduct);}
1262
operator ^(OperatorOnUnknown & opu,const real_t & val)1263 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const real_t& val)
1264 {return updateRight(opu, val, _crossProduct);}
1265
operator %(OperatorOnUnknown & opu,const real_t & val)1266 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const real_t& val)
1267 {return updateRight(opu, val, _contractedProduct);}
1268
operator *(const real_t & val,OperatorOnUnknown & opu)1269 OperatorOnUnknown& operator*(const real_t& val, OperatorOnUnknown& opu)
1270 {return updateLeft(opu, val, _product);}
1271
operator |(const real_t & val,OperatorOnUnknown & opu)1272 OperatorOnUnknown& operator|(const real_t& val, OperatorOnUnknown& opu)
1273 {return updateLeft(opu, val, _innerProduct);}
1274
operator ^(const real_t & val,OperatorOnUnknown & opu)1275 OperatorOnUnknown& operator^(const real_t& val, OperatorOnUnknown& opu)
1276 {return updateLeft(opu, val, _crossProduct);}
1277
operator %(const real_t & val,OperatorOnUnknown & opu)1278 OperatorOnUnknown& operator%(const real_t& val, OperatorOnUnknown& opu)
1279 {return updateLeft(opu, val, _contractedProduct);}
1280
1281 // left and right operation with complex scalar : op(u) aop c or c aop op(u)
operator *(OperatorOnUnknown & opu,const complex_t & val)1282 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const complex_t& val)
1283 {return updateRight(opu, val, _product);}
1284
operator |(OperatorOnUnknown & opu,const complex_t & val)1285 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const complex_t& val)
1286 {return updateRight(opu, val, _innerProduct);}
1287
operator ^(OperatorOnUnknown & opu,const complex_t & val)1288 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const complex_t& val)
1289 {return updateRight(opu, val, _crossProduct);}
1290
operator %(OperatorOnUnknown & opu,const complex_t & val)1291 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const complex_t& val)
1292 {return updateRight(opu, val, _contractedProduct);}
1293
operator *(const complex_t & val,OperatorOnUnknown & opu)1294 OperatorOnUnknown& operator*(const complex_t& val, OperatorOnUnknown& opu)
1295 {return updateLeft(opu, val, _product);}
1296
operator |(const complex_t & val,OperatorOnUnknown & opu)1297 OperatorOnUnknown& operator|(const complex_t& val, OperatorOnUnknown& opu)
1298 {return updateLeft(opu, val, _innerProduct);}
1299
operator ^(const complex_t & val,OperatorOnUnknown & opu)1300 OperatorOnUnknown& operator^(const complex_t& val, OperatorOnUnknown& opu)
1301 {return updateLeft(opu, val, _crossProduct);}
1302
operator %(const complex_t & val,OperatorOnUnknown & opu)1303 OperatorOnUnknown& operator%(const complex_t& val, OperatorOnUnknown& opu)
1304 {return updateLeft(opu, val, _contractedProduct);}
1305
1306 //==========================================================================================
1307 // operators involving OperatorOnUnknown and Value (update current OperatorOnUnknown)
1308 //==========================================================================================
1309 // update left and right operation with Value : op(u) aop Value or Value aop op(u)
updateRight(OperatorOnUnknown & opu,const Value & val,AlgebraicOperator o)1310 OperatorOnUnknown& updateRight(OperatorOnUnknown& opu, const Value& val, AlgebraicOperator o)
1311 {
1312 if(opu.rightOperand() != 0) { error("operator_sideisdef", words("right")); }
1313 //opu.updateReturnedType(o, val.valueType(), val.strucType(), val.dims(), false);
1314 opu.rightOperand() = new Operand(&val, o);
1315 if(opu.leftOperand() == 0) { opu.leftPriority() = false; }
1316 opu.setStructure();
1317 return opu;
1318 }
updateLeft(OperatorOnUnknown & opu,const Value & val,AlgebraicOperator o)1319 OperatorOnUnknown& updateLeft(OperatorOnUnknown& opu, const Value& val, AlgebraicOperator o)
1320 {
1321 if(opu.leftOperand() != 0) { error("operator_sideisdef", words("left")); }
1322 //opu.updateReturnedType(o, val.valueType(), val.strucType(), val.dims(), true);
1323 opu.leftOperand() = new Operand(&val, o);
1324 if(opu.rightOperand() == 0) { opu.leftPriority() = true; }
1325 opu.setStructure();
1326 return opu;
1327 }
operator *(OperatorOnUnknown & opu,const Value & val)1328 OperatorOnUnknown& operator*(OperatorOnUnknown& opu, const Value& val) // d(u) * val
1329 {
1330 return updateRight(opu, val, _product);
1331 }
operator |(OperatorOnUnknown & opu,const Value & val)1332 OperatorOnUnknown& operator|(OperatorOnUnknown& opu, const Value& val) // d(u) | val
1333 {
1334 return updateRight(opu, val, _innerProduct);
1335 }
operator ^(OperatorOnUnknown & opu,const Value & val)1336 OperatorOnUnknown& operator^(OperatorOnUnknown& opu, const Value& val) // d(u) ^ val
1337 {
1338 return updateRight(opu, val, _crossProduct);
1339 }
operator %(OperatorOnUnknown & opu,const Value & val)1340 OperatorOnUnknown& operator%(OperatorOnUnknown& opu, const Value& val) // d(u) % val
1341 {
1342 return updateRight(opu, val, _contractedProduct);
1343 }
operator *(const Value & val,OperatorOnUnknown & opu)1344 OperatorOnUnknown& operator*(const Value& val, OperatorOnUnknown& opu) // val * d(u)
1345 {
1346 return updateLeft(opu, val, _product);
1347 }
operator |(const Value & val,OperatorOnUnknown & opu)1348 OperatorOnUnknown& operator|(const Value& val, OperatorOnUnknown& opu) // val | d(u)
1349 {
1350 return updateLeft(opu, val, _innerProduct);
1351 }
operator ^(const Value & val,OperatorOnUnknown & opu)1352 OperatorOnUnknown& operator^(const Value& val, OperatorOnUnknown& opu) // val ^ d(u)
1353 {
1354 return updateLeft(opu, val, _crossProduct);
1355 }
operator %(const Value & val,OperatorOnUnknown & opu)1356 OperatorOnUnknown& operator%(const Value& val, OperatorOnUnknown& opu) // val % d(u)
1357 {
1358 return updateLeft(opu, val, _contractedProduct);
1359 }
1360
1361 //==========================================================================================
1362 // various facilities
1363 //==========================================================================================
1364
1365 // check is operation between two operatoronunknown is consistant in terms of structure
1366 // the result have to be a scalar
checkConsistancy(const OperatorOnUnknown & opu,AlgebraicOperator aop,const OperatorOnUnknown & opv)1367 bool checkConsistancy(const OperatorOnUnknown& opu, AlgebraicOperator aop, const OperatorOnUnknown& opv)
1368 {
1369 if(opu.strucType() == _scalar && opv.strucType() == _scalar) { return true; }
1370 switch(aop)
1371 {
1372 case _product : // only scalar * scalar
1373 if(opu.strucType() == _scalar && opv.strucType() == _scalar) { return true; }
1374 break;
1375 case _innerProduct : // only vector | vector
1376 if(opu.strucType() == _vector && opv.strucType() == _vector) { return true; }
1377 break;
1378 case _contractedProduct : // only matrix % matrix
1379 if(opu.strucType() == _matrix && opv.strucType() == _matrix) { return true; }
1380 break;
1381 case _crossProduct : // never
1382 break;
1383 }
1384 return false;
1385 }
1386
1387
1388 // compare OperatorOnUnknown (same unknowns, same diff operators, same functions ...)
operator ==(const OperatorOnUnknown & opu,const OperatorOnUnknown & opv)1389 bool operator==(const OperatorOnUnknown& opu, const OperatorOnUnknown& opv)
1390 {
1391 if(opu.unknown()!=opv.unknown()) return false;
1392 if(opu.conjugateUnknown()!=opv.conjugateUnknown()) return false;
1393 if(opu.difOp()!=opv.difOp()) return false;
1394 if(opu.leftOperand()!=0 && opv.leftOperand()!=0)
1395 if(*opu.leftOperand()!=*opv.leftOperand()) return false;
1396 if(opu.rightOperand()!=0 && opv.rightOperand()!=0)
1397 if(*opu.rightOperand()!=*opv.rightOperand()) return false;
1398 if(opu.coefs()!=opu.coefs()) return false;
1399 return true;
1400 }
1401
operator !=(const OperatorOnUnknown & opu,const OperatorOnUnknown & opv)1402 bool operator!=(const OperatorOnUnknown& opu, const OperatorOnUnknown& opv)
1403 {
1404 return !(opu==opv);
1405 }
1406
1407 //==========================================================================================
1408 // printing facilities
1409 //==========================================================================================
print(std::ostream & os) const1410 void OperatorOnUnknown::print(std::ostream& os) const
1411 {
1412 if(theVerboseLevel > 0)
1413 {
1414 os << words("operator") << " " << difOp_p->name() << " "<<u_p->name() ;
1415 if(conjugateUnknown_) { os << "(" << words("conjugate") << ")"; }
1416 os << " " << words("returns a") << " " << words("value", type_) << " " << words("structure", struct_);
1417 if(struct_!=_scalar) os << " ("<<dimsRes_.first<<" x "<<dimsRes_.second<<")";
1418 if(coefs_.size()>0) os<<" coefficients : "<<coefs_;
1419 if(leftOperand_p != 0) { os << "\n " << words("left operand") << " : " << *leftOperand_p; }
1420 if(rightOperand_p != 0) { os << "\n " << words("right operand") << " : " << *rightOperand_p; }
1421 if(leftOperand_p != 0 && rightOperand_p != 0)
1422 {
1423 if(leftPriority_) { os << " (" << words("left priority") << ")"; }
1424 else { os << " (" << words("right priority") << ")"; }
1425 }
1426 os << "\n";
1427 }
1428 }
1429
asString() const1430 string_t OperatorOnUnknown::asString() const
1431 {
1432 string_t s="";
1433 if(leftOperand_p != 0) s += leftOperand_p->asString();
1434 if(difOp_p->type()!=_id) s += difOp_p->name()+"("+u_p->name()+")";
1435 else s += u_p->name();
1436 if(rightOperand_p != 0) s += rightOperand_p->asString();
1437 return s;
1438 }
1439
printsymbolic(std::ostream & os) const1440 void OperatorOnUnknown::printsymbolic(std::ostream& os) const
1441 {
1442 if(leftOperand_p != 0) { leftOperand_p->printsymbolic(os);}
1443 if(difOp_p->type()!=_id)
1444 os<<difOp_p->name()<<"("<<u_p->name()<<")";
1445 else
1446 os<<u_p->name();
1447 if(rightOperand_p != 0) { rightOperand_p->printsymbolic(os);}
1448 }
1449
operator <<(std::ostream & os,const OperatorOnUnknown & opu)1450 std::ostream& operator<<(std::ostream& os, const OperatorOnUnknown& opu)
1451 {
1452 opu.print(os);
1453 return os;
1454 }
1455
1456 } // end of namespace xlifepp
1457
1458