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 termUtils.cpp
19 \author N. Kielbasiewicz
20 \since 27 feb 2013
21 \date 1 mar 2013
22
23 \brief Implementation of term utility functions
24 */
25
26 #include "termUtils.hpp"
27 #include "../TermMatrix.hpp"
28 #include "../KernelOperatorOnTermVector.hpp"
29 #include "../../form/BilinearForm.hpp"
30 #include "SauterSchwabIM.hpp"
31 #include "DuffyIM.hpp"
32 #include <vector>
33 #include <sstream>
34
35 namespace xlifepp
36 {
37
38 //! special cast function used by FEComputation's
toRealComplex(const Vector<real_t> & v,const complex_t & c)39 real_t toRealComplex(const Vector<real_t>& v,const complex_t& c) {return c.real();}
toRealComplex(const Vector<complex_t> &,const complex_t & c)40 complex_t toRealComplex(const Vector<complex_t>&,const complex_t& c) {return c;}
toRealComplex(const Vector<Vector<real_t>> & v,const complex_t & c)41 real_t toRealComplex(const Vector<Vector<real_t> >& v,const complex_t& c) {return c.real();}
toRealComplex(const Vector<Vector<complex_t>> & v,const complex_t & c)42 complex_t toRealComplex(const Vector<Vector<complex_t> >& v,const complex_t& c) {return c;}
43
44 // build a storage from a pair of Space
45 // for FE space it is assumed that only connected dofs in contiguous elements give non zero values
46 // FE spaces must live on same domain
buildStorage(const Space & rs,const Space & cs,StorageType st,AccessType at)47 MatrixStorage* buildStorage(const Space& rs, const Space& cs, StorageType st, AccessType at)
48 {
49 trace_p->push("buildStorage");
50
51 // the storage id is built from the adresses of Space objects
52 std::stringstream ss;
53 ss << &rs << "-" << &cs;
54 MatrixStorage* res = findMatrixStorage(ss.str(), st, at);
55
56 if(res != 0)
57 {
58 // the storage already exists, then we return a pointer to
59 trace_p->pop();
60 return res;
61 }
62
63 // we have to build a new storage
64 //detect special case : one domain is the extension of the other
65 const MeshDomain* cdom = cs.domain()->meshDomain(), *rdom = rs.domain()->meshDomain();
66 string_t extension="";
67 if(cdom!=0 && rdom!=0 && cdom!=rdom)
68 {
69 if(cdom->extensionof_p == rdom) extension="c=ext_r";
70 if(rdom->extensionof_p == cdom) extension="r=ext_c";
71 }
72 std::vector<number_t>::const_iterator cit_r, cit_c;
73
74 //build col index
75 number_t dimr=rs.nbDofs(), dimc=cs.nbDofs();
76 std::vector<std::set<number_t> > colset; // temporary vector of col index
77 if(st!=_dense)
78 {
79 colset.resize(dimr);
80 if(rs.typeOfSubSpace() == _spSpace || cs.typeOfSubSpace() == _spSpace || (cs.domain()!=rs.domain() && extension==""))
81 {
82 std::set<number_t> colInd;
83 std::set<number_t>::iterator jt=colInd.begin();
84 for(number_t j = 0; j < dimc; j++) jt=colInd.insert(jt, j + 1);
85 std::vector<std::set<number_t> >::iterator it=colset.begin();
86 for(number_t i = 0; i < dimr; i++, it++) { *it = colInd; }
87 }
88
89 else if(extension=="c=ext_r") // case of extension of row domain
90 {
91
92 std::set<GeomElement*> elts_ext(cdom->geomElements.begin(),cdom->geomElements.end()); //set of column geom elements
93 number_t n = rs.nbOfElements();
94 for(number_t k = 0; k < n; k++) //loop on side domain elements
95 {
96 std::vector<number_t> red = rs.elementDofs(k); // row dof ranks for element k
97 GeomElement* gelt = rs.element_p(k)->geomElt_p; // geometric element supporting side element
98 std::vector<GeoNumPair>& parents=gelt->parentSides(); // parent of side element
99 std::vector<GeoNumPair>::iterator itp=parents.begin();
100 for(; itp!=parents.end(); ++itp) //loop on parents
101 if(elts_ext.find(itp->first)!=elts_ext.end()) //parent elts in extended domain
102 {
103 std::vector<number_t> ced = cs.elementDofs(cs.numElement(itp->first)); // col dof ranks for element k ### A VERIFIER
104 for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
105 for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
106 colset[*cit_r - 1].insert(*cit_c);
107 }
108 }
109 }
110 else if(extension=="r=ext_c") // case of a extension of col domain
111 {
112 std::set<GeomElement*> elts_ext(rdom->geomElements.begin(),rdom->geomElements.end());//set of row geom elements
113 number_t n = cs.nbOfElements();
114 for(number_t k = 0; k < n; k++) //loop on side domain elements
115 {
116 std::vector<number_t> ced = cs.elementDofs(k); // col dof ranks for element k
117 GeomElement* gelt = cs.element_p(k)->geomElt_p; // geometric element supporting side element
118 std::vector<GeoNumPair>& parents=gelt->parentSides(); // parent of side element
119 std::vector<GeoNumPair>::iterator itp=parents.begin();
120 for(; itp!=parents.end(); ++itp) //loop on parents
121 if(elts_ext.find(itp->first)!=elts_ext.end()) //parent elts in extended domain
122 {
123 std::vector<number_t> red = rs.elementDofs(rs.numElement(itp->first)); // row dof ranks for element k ### A VERIFIER
124 for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
125 for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
126 colset[*cit_r - 1].insert(*cit_c);
127 }
128 }
129 }
130 else // both spaces are standard finite element spaces
131 {
132 number_t n = rs.nbOfElements();
133 // build list of col index for each rows
134 if(at!=_sym)
135 {
136 for(number_t k = 0; k < n; k++)
137 {
138 // dof ranks for element k
139 std::vector<number_t> red = rs.elementDofs(k); // row dof ranks for element k
140 GeomElement* gelt = rs.element_p(k)->geomElt_p; // geometric element supporting element
141 number_t kc=k;
142 if(cs.element_p(k)->geomElt_p!=gelt) kc=cs.numElement(gelt); //not the same support, find the right one
143 std::vector<number_t> ced = cs.elementDofs(kc); // col dof ranks for element k
144 for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
145 for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
146 colset[*cit_r - 1].insert(*cit_c);
147 }
148 }
149 else //symmetric storage : keep only strict lower part indices
150 {
151 for(number_t k = 0; k < n; k++)
152 {
153 std::vector<number_t> red = rs.elementDofs(k); // row dof ranks for element k
154 GeomElement* gelt = rs.element_p(k)->geomElt_p; // geometric element supporting element
155 number_t kc=k;
156 if(cs.element_p(k)->geomElt_p!=gelt) kc=cs.numElement(gelt); //not the same support, find the right one
157 std::vector<number_t> ced = cs.elementDofs(kc); // col dof ranks for element k
158 for(cit_r = red.begin(); cit_r != red.end(); cit_r++)
159 for(cit_c = ced.begin(); cit_c != ced.end(); cit_c++)
160 if(*cit_r > *cit_c) colset[*cit_r - 1].insert(*cit_c);
161 }
162 }
163 }
164 }
165 // construct the storage with the column indices by rows
166 switch(st)
167 {
168 case _cs :
169 switch(at)
170 {
171 case _sym : trace_p->pop(); return new SymCsStorage(dimr, colset, _lower, ss.str());
172 case _row : trace_p->pop(); return new RowCsStorage(dimr, dimc, colset, ss.str());
173 case _col : trace_p->pop(); return new ColCsStorage(dimr, dimc, colset, ss.str());
174 case _dual : trace_p->pop(); return new DualCsStorage(dimr, dimc, colset, ss.str());
175 default : error("storage_bad_access", words("access type",at), words("storage type",_cs));
176 }
177 break;
178 case _skyline :
179 switch(at)
180 {
181 case _sym : trace_p->pop(); return new SymSkylineStorage(dimr, colset, ss.str());
182 case _dual : trace_p->pop(); return new DualSkylineStorage(dimr, dimc, colset, ss.str());
183 default : error("storage_bad_access", words("access type",at), words("storage type",_skyline));
184 }
185 break;
186 case _dense :
187 switch(at)
188 {
189 case _sym : trace_p->pop(); return new SymDenseStorage(dimr, ss.str());
190 case _row : trace_p->pop(); return new RowDenseStorage(dimr, dimc, ss.str());
191 case _col : trace_p->pop(); return new ColDenseStorage(dimr, dimc, ss.str());
192 case _dual : trace_p->pop(); return new DualDenseStorage(dimr, dimc, ss.str());
193 default : error("storage_bad_access", words("access type",at), words("storage type",_dense));
194 }
195 break;
196 default : error("storage_not_implemented","buildStorage",words("storage type",st));
197 }
198 trace_p->pop();
199 return 0;
200 }
201
202 /*! modify MatrixStorage stm to include dense submatrix given by its row and column indices (rows and cols)
203 the new storage is of type (st,at)
204 rows (resp. cols) has to be a subset of row indices (resp. column indices) of the stm storage
205 if overwrite=true the current storage is overwritten and returned
206 else a new MatrixStorage is returned, the stm storage being not cleared!
207 */
208
updateStorage(MatrixStorage & stm,const std::vector<number_t> rows,const std::vector<number_t> cols,StorageType st,AccessType at,bool overwrite)209 MatrixStorage* updateStorage(MatrixStorage& stm, const std::vector<number_t> rows, const std::vector<number_t> cols,
210 StorageType st, AccessType at, bool overwrite)
211 {
212 trace_p->push("updateStorage(MatrixStorage, rows, cols,...)");
213 MatrixStorage* newstm=&stm;
214 if(!overwrite) newstm= stm.clone();
215 if(stm.storageType()==st && stm.accessType()==at) //target storage is same as stm storage
216 {
217 switch(st)
218 {
219 case _cs :
220 case _skyline : newstm->addSubMatrixIndices(rows,cols);
221 default : //dense in dense, nothing to do
222 trace_p->pop();
223 return newstm;
224 }
225 }
226
227 error("storage_not_handled", words("storage type",stm.storageType()),words("access type",stm.accessType()));
228 trace_p->pop();
229 return newstm;
230 }
231
232 // build a storage from type and dimension, no allocation of pointers (void matrix)
buildStorage(StorageType st,AccessType at,number_t dimr,number_t dimc)233 MatrixStorage* buildStorage(StorageType st, AccessType at, number_t dimr, number_t dimc)
234 {
235 switch(st)
236 {
237 case _cs :
238 switch(at)
239 {
240 case _sym : return new SymCsStorage(dimr);
241 case _row : return new RowCsStorage(dimr, dimc);
242 case _col : return new ColCsStorage(dimr, dimc);
243 case _dual : return new DualCsStorage(dimr, dimc);
244 default : error("storage_bad_access", words("access type",at), words("storage type",_cs));
245 }
246 break;
247 case _skyline :
248 switch(at)
249 {
250 case _sym : return new SymSkylineStorage(dimr);
251 case _dual : return new DualSkylineStorage(dimr, dimc);
252 default : error("storage_bad_access", words("access type",at), words("storage type",_skyline));
253 }
254 break;
255 case _dense :
256 switch(at)
257 {
258 case _sym : return new SymDenseStorage(dimr);
259 case _row : return new RowDenseStorage(dimr, dimc);
260 case _col : return new ColDenseStorage(dimr, dimc);
261 case _dual : return new DualDenseStorage(dimr, dimc);
262 default : error("storage_bad_access", words("access type",at), words("storage type",_dense));
263 }
264 break;
265 default : error("storage_not_implemented","buildStorage",words("storage type",st));
266 }
267 return 0;
268 }
269
270 // build a storage from type, dimension and column indices (vector of vectors)
buildStorage(StorageType st,AccessType at,number_t dimr,number_t dimc,const std::vector<std::vector<number_t>> & indices)271 MatrixStorage* buildStorage(StorageType st, AccessType at, number_t dimr, number_t dimc, const std::vector<std::vector<number_t> >& indices)
272 {
273 trace_p->push("buildStorage");
274 string_t str="";
275 switch(st)
276 {
277 case _cs :
278 switch(at)
279 {
280 case _sym : trace_p->pop(); return new SymCsStorage(dimr, indices, _all, str);
281 case _row : trace_p->pop(); return new RowCsStorage(dimr, dimc, indices, str);
282 case _col :
283 {
284 trace_p->pop();
285 std::vector<std::vector<number_t> > rowindices(dimc);
286 std::vector<std::vector<number_t> >::const_iterator itr=indices.begin();
287 std::vector<number_t>::const_iterator itc;
288 number_t r=1;
289 for(; itr!=indices.end(); itr++, r++)
290 for(itc=itr->begin(); itc!= itr->end(); itc++) rowindices[*itc - 1].push_back(r);
291 return new ColCsStorage(dimr, dimc, rowindices, str);
292 }
293 case _dual : trace_p->pop(); return new DualCsStorage(dimr, dimc, indices, str);
294 default : error("storage_bad_access", words("access type",at), words("storage type",_cs));
295 }
296 break;
297 case _skyline :
298 switch(at)
299 {
300 case _sym : trace_p->pop(); return new SymSkylineStorage(dimr, indices, str);
301 case _dual : trace_p->pop(); return new DualSkylineStorage(dimr, dimc, indices, str);
302 default : error("storage_bad_access", words("access type",at), words("storage type",_skyline));
303 }
304 break;
305 case _dense :
306 switch(at)
307 {
308 case _sym : trace_p->pop(); return new SymDenseStorage(dimr, str);
309 case _row : trace_p->pop(); return new RowDenseStorage(dimr, dimc, str);
310 case _col : trace_p->pop(); return new ColDenseStorage(dimr, dimc, str);
311 case _dual : trace_p->pop(); return new DualDenseStorage(dimr, dimc, str);
312 default : error("storage_bad_access", words("access type",at), words("storage type",_dense));
313 }
314 break;
315 default : error("storage_not_implemented","buildStorage",words("storage type",st));
316 }
317
318 trace_p->pop();
319 return 0;
320 }
321
322 // build a storage from type, dimension and column indices (vector of sets)
buildStorage(StorageType st,AccessType at,number_t dimr,number_t dimc,const std::vector<std::set<number_t>> & colIndices)323 MatrixStorage* buildStorage(StorageType st, AccessType at, number_t dimr, number_t dimc, const std::vector<std::set<number_t> >& colIndices)
324 {
325 // convert set to vector
326 std::vector<std::vector<number_t> > indices(colIndices.size());
327 std::vector<std::set<number_t> >::const_iterator itc;
328 std::vector<std::vector<number_t> >::iterator iti = indices.begin();
329 for(itc = colIndices.begin(); itc != colIndices.end(); ++itc, ++iti)
330 *iti = std::vector<number_t>(itc->begin(), itc->end());
331 return buildStorage(st,at,dimr,dimc,indices);
332 }
333
334 /*! compute interpolated normals on Lagrange Dofs of space sp
335 solve the projection problem Mk*nk = M0*n0 where
336 Mk is the mass matrix associated to sp interpolation
337 M0 is the "hybrid" mass matrix between sp interpolation and P0 interpolation
338 the normals are returned as a vector of Vector
339
340 NOTE : this method is well adapted for planar element (same normal every where in element)
341 in case of curved geometric element, the method may be improved by using Gauss Lobatto interpolation of higher degree
342 */
interpolatedNormals(Space & sp,std::vector<Vector<real_t>> & ns)343 void interpolatedNormals(Space& sp, std::vector<Vector<real_t> >& ns)
344 {
345 if(!sp.isFE())
346 {
347 where("interpolatedNormals(Space, vector<Vector<Real> >");
348 error("not_fe_space_type", sp.name());
349 }
350 const Interpolation* intp=sp.rootSpace()->interpolation();
351 if(intp->type!=_Lagrange) error("lagrange_fe_space_only",sp.name());
352
353 trace_p->push("interpolatedNormals");
354
355 if(intp->numtype > 0) //non P0 case
356 {
357
358 const GeomDomain& dom=*sp.domain();
359 dimen_t d=dom.mesh()->spaceDim();
360 Unknown n(*sp.rootSpace(),"n",d);
361 TestFunction t(n,"t");
362
363 //define matrix
364 TermMatrix Mk(intg(dom,n|t),"normalPkM");
365 Space V0(dom,P0,"V0",false);
366 Unknown n0(V0,"n0",d);
367 TermMatrix M0(intg(dom,n0|t),"normalPkP0M");
368
369 //fill normal P0 in a TermVector
370 TermVector N0(n0, dom, Vector<real_t>(d,0.),"normalsP0");
371 Vector<Vector<real_t> >::iterator it=N0.begin()->second->entries()->beginrv();
372 for(number_t k=0; k<sp.nbOfElements(); k++, it++)
373 {
374 const Element* elt=sp.element_p(k);
375 MeshElement* melt = elt->geomElt_p->meshElement();
376 if(melt==0) melt=elt->geomElt_p->buildSideMeshElement();
377 GeomMapData gmap(melt);
378 gmap.computeJacobianMatrix(elt->refElt_p->geomRefElem_p->centroid());
379 gmap.computeOutwardNormal();
380 *it = gmap.normalVector;
381 }
382
383 //solve projection problem
384 TermVector Nk=directSolve(Mk,M0*N0);
385
386 //return normals
387 ns= *Nk.begin()->second->entries()->rvEntries_p;
388 std::vector<Vector<real_t> >::iterator itn=ns.begin();
389 for(; itn!=ns.end(); itn++) *itn/=itn->norm2();
390 }
391
392 else //P0 interpolation
393 {
394 verboseLevel(10);
395 ns.resize(sp.nbOfElements());
396 std::vector<Vector<real_t> >::iterator it=ns.begin();
397 for(number_t k=0; k<sp.nbOfElements(); k++, it++)
398 {
399 const Element* elt=sp.element_p(k);
400 MeshElement* melt = elt->geomElt_p->meshElement();
401 if(melt==0) melt=elt->geomElt_p->buildSideMeshElement();
402 GeomMapData gmap(melt);
403 gmap.computeJacobianMatrix(elt->refElt_p->geomRefElem_p->centroid());
404 gmap.computeOutwardNormal();
405 *it = gmap.normalVector;
406 }
407 }
408 trace_p->pop();
409 return;
410 }
411
412 //compute normals of a side domain using interpolation given by an unknown
normalsOn(GeomDomain & dom,const Unknown & u)413 TermVector normalsOn(GeomDomain& dom, const Unknown& u)
414 {
415 if (!dom.isSideDomain() && ((dom.dim()+1)!=dom.spaceDim()))
416 {
417 where("normalsOn(Domain, Unknown)");
418 error("domain_not_side", dom.name());
419 }
420 const Space* spu = u.space();
421 const Interpolation* itp0 = findInterpolation(Lagrange, _standard, 0, L2);
422 const Space* sp0=0; bool sub = false;
423 if (spu->interpolation()==itp0)
424 {
425 if (spu->domain() == &dom) sp0 = spu; // u space is a P0 space
426 else if (spu->domain()->include(dom)) {sp0 = spu; sub=true;}
427 }
428 if (sp0==0) // find or create space P0
429 {
430 std::vector<Space*>::iterator it;
431 for (it = Space::theSpaces.begin(); it != Space::theSpaces.end() && (sp0==0 || sub==false); ++it)
432 {
433 const Space* spit = (*it)->space();
434 if (spit != 0)
435 {
436 const Interpolation* itp = spit->interpolation();
437 if (itp == itp0)
438 {
439 if (spit->domain()==&dom) sp0=spit;
440 else if (spit->domain()->include(dom)) {sp0=spit; sub=true;}
441 }
442 }
443 }
444 }
445 if (sp0==0) sp0 = new Space(dom, P0, "P0"); // no P0 space found
446 else if (sub) sp0 = Space::findSubSpace(&dom, sp0); // P0 parent space found
447 dimen_t d = dom.spaceDim();
448 Unknown n0(*const_cast<Space*>(sp0), "n_"+dom.name(),d);
449
450 // fill normal P0 in a TermVector
451 TermVector N0(n0, dom, Vector<real_t>(d,0.),"normalsP0");
452 Vector<Vector<real_t> >::iterator it=N0.begin()->second->entries()->beginrv();
453 for (number_t k=0; k<sp0->nbOfElements(); k++, it++)
454 {
455 if (!dom.meshDomain()->normalComputed)
456 {
457 const Element* elt=sp0->element_p(k);
458 MeshElement* melt = elt->geomElt_p->meshElement();
459 if (melt==0) melt=elt->geomElt_p->buildSideMeshElement();
460 GeomMapData gmap(melt);
461 gmap.computeJacobianMatrix(elt->refElt_p->geomRefElem_p->centroid());
462 gmap.computeOutwardNormal();
463 *it = gmap.normalVector;
464 }
465 else *it = sp0->element_p(k)->geomElt_p->meshElement()->geomMapData_p->normalVector;
466
467 }
468
469 if (spu->interpolation()==itp0) return N0;
470
471 //project normal on space spu
472 Unknown* uc = const_cast<Unknown*>(&u);
473 //if(spu->dimFun()==1 && u.nbOfComponents()==1) uc=new Unknown(*const_cast<Space*>(spu),"uc_"+tostring(spu),d); // permissive, user pass a scalar unknown
474 if (uc->nbOfComponents()!=d && spu->dimFun()!=d) error("unknown_inconsistent_size", uc->name(), d, uc->nbOfComponents());
475 TestFunction t(*uc,"t_u");
476 TermMatrix Mk(intg(dom,(*uc)|t),"u_t_mass");
477 TermMatrix M0(intg(dom,n0|t),"n0_t_mass");
478 TermVector Nk=directSolve(Mk,M0*N0);
479 return Nk;
480 }
481
482 /*! precompute geometry data and shape values on any element for any quadrature involved
483 restricted to one order geometric element, take large place but computation is then faster in case of IE!!!
484
485 dom : integration domain
486 subsp : subspace providing the list of elements
487 quads : list of quadratures
488 nbc : number of components of the unknown
489 der : compute shape derivatives if true
490 nor : compute normal if true
491 extendedDomain : true when dealing with an extended domain
492 */
preComputationIE(const GeomDomain & dom,const Space * subsp,std::set<Quadrature * > & quads,dimen_t nbc,bool der,bool nor,bool extendedDomain)493 void preComputationIE(const GeomDomain& dom, const Space* subsp, std::set<Quadrature*>& quads, dimen_t nbc, bool der, bool nor,
494 bool extendedDomain)
495 {
496 der=true; //force computation of derivatives to avoid misfit when recall
497
498 //shapevalues for all single quadrature (non parallel)
499 for(number_t k = 0; k < subsp->nbOfElements(); k++) // data related to element
500 {
501 const Element* elt = subsp->element_p(k);
502 RefElement* relt = elt->refElt_p; //reference element related to element
503 GeomElement* gelt = elt->geomElt_p;
504 const GeomRefElement* grelt = gelt->geomRefElement();
505 number_t sw = relt->nbDofs()*nbc;
506 if(extendedDomain) //build geometric data for extended element
507 {
508 if(gelt->meshElement()==0) gelt->buildSideMeshElement();
509 MeshElement* melt = gelt->meshElement();
510 GeomMapData* mapdata = melt->geomMapData_p;
511 if(mapdata==0) mapdata=new GeomMapData(melt);
512 mapdata->computeJacobianMatrix(std::vector<real_t>(gelt->elementDim(),0.));
513 mapdata->invertJacobianMatrix();
514 melt->geomMapData_p = mapdata;
515 }
516 //precompute reference shape values
517 std::set<Quadrature*>::iterator itq=quads.begin();
518 for(; itq!=quads.end(); ++itq)
519 {
520 std::map<Quadrature*,std::vector<ShapeValues> >::iterator itmq=relt->qshvs_.find(*itq);
521 if(itmq==relt->qshvs_.end() || itmq->second.begin()->w.size()!=sw) //shapevalues re-computed
522 {
523 dimen_t dq = (*itq)->dim(); //dimension of quadrature points
524 number_t nq = (*itq)->numberOfPoints(); //number of quadrature points
525 if(!extendedDomain) //no extension
526 {
527 std::vector<real_t>::const_iterator itp = (*itq)->point();
528 relt->qshvs_[*itq]=std::vector<ShapeValues>(nq); //to store all shapevalues at quadrature points
529 std::vector<ShapeValues>::iterator itshv = relt->qshvs_[*itq].begin();
530 for(number_t i = 0; i < nq; i++, itp += dq, itshv++)
531 {
532 relt->computeShapeValues(itp, der);
533 *itshv = relt->shapeValues;
534 if(nbc>1) itshv->extendToVector(nbc); //extend scalar shape functions to nbc vector shape functions
535 }
536 }
537 else // in case of extension build shapevalues for each side of ref element stored sequentially
538 {
539 number_t nbs = relt->geomRefElem_p->nbSides();
540 relt->qshvs_[*itq]=std::vector<ShapeValues>(nbs*nq); //to store all shapevalues at quadrature points
541 std::vector<ShapeValues>::iterator itshv = relt->qshvs_[*itq].begin();
542 for(number_t s=1; s<=nbs; ++s) //loop on side
543 {
544 std::vector<real_t>::const_iterator itp = (*itq)->point();
545 for(number_t i = 0; i < nq; i++, itp += dq, itshv++)
546 {
547 relt->computeShapeValues(grelt->sideToElt(s,itp).begin(), der);
548 *itshv = relt->shapeValues;
549 if(nbc>1) itshv->extendToVector(nbc); //extend scalar shape functions to nbc vector shape functions
550 }
551 }
552 }
553 }
554 }
555 }
556 //precompute signDofs of element if required
557 #ifdef XLIFEPP_WITH_OMP
558 #pragma omp parallel for
559 #endif // XLIFEPP_WITH_OMP
560 for(number_t k = 0; k < subsp->nbOfElements(); k++) // data related to element
561 {
562 const Element* elt = subsp->element_p(k);
563 if(elt->refElt_p->dofCompatibility == _signDofCompatibility) elt->setDofSigns();
564 }
565
566 MeshDomain* mdom=const_cast<MeshDomain*>(dom.meshDomain());
567 //compute map of quadrature points (parallel)
568 #ifdef XLIFEPP_WITH_OMP
569 #pragma omp parallel for
570 #endif // XLIFEPP_WITH_OMP
571 for(number_t k = 1; k <= mdom->numberOfElements(); k++)
572 {
573 GeomElement* gelt =mdom->element(k);
574 //geometry information
575 if(gelt->meshElement()==0) gelt->buildSideMeshElement();
576 MeshElement* melt = gelt->meshElement();
577 GeomMapData* mapdata=melt->geomMapData_p;
578 if(mapdata==0) mapdata=new GeomMapData(melt);
579 mapdata->computeJacobianMatrix(std::vector<real_t>(gelt->elementDim(),0.));
580 mapdata->invertJacobianMatrix();
581 mapdata->computeDifferentialElement();
582 melt->geomMapData_p = mapdata;
583 if(nor) mapdata->computeOutwardNormal();
584 //physical points
585 std::set<Quadrature*>::iterator itq=quads.begin();
586 for(; itq!=quads.end(); ++itq)
587 {
588 dimen_t dq = (*itq)->dim(); //dimension of quadrature points
589 number_t nq = (*itq)->numberOfPoints(); //number of quadrature points
590 std::map<Quadrature*,std::vector<Point> >::iterator itmqphy= mapdata->phyPoints.find(*itq);
591 if(itmqphy==mapdata->phyPoints.end()) //create physical points
592 {
593 mapdata->phyPoints[*itq] = std::vector<Point>(nq);
594 std::vector<Point>::iterator itphy=mapdata->phyPoints[*itq].begin();
595 std::vector<real_t>::const_iterator itp = (*itq)->point();
596 for(number_t i = 0; i < nq; i++, itp += dq, ++itphy) *itphy = mapdata->geomMap(itp);
597 }
598 }
599 }
600 }
601
602 /*! precompute geometry data and shape values on any element for any quadrature involved
603 restricted to one order geometric element, take large place but computation is then faster in case of IE!!!
604
605 subsp : subspace providing the list of elements
606 quads : list of quadratures
607 nbc : number of components of the unknown
608 der : compute shape derivatives if true
609 nor : compute normal if true
610 */
preComputationIE_old(const Space * subsp,std::set<Quadrature * > & quads,dimen_t nbc,bool der,bool nor)611 void preComputationIE_old(const Space* subsp, std::set<Quadrature*>& quads, dimen_t nbc, bool der, bool nor)
612 {
613 der=true; //force computation of derivatives to avoid misfit when recall
614
615 for(number_t k = 0; k < subsp->nbOfElements(); k++) // data related to u
616 {
617 const Element* elt = subsp->element_p(k);
618 GeomElement* gelt = elt->geomElt_p;
619
620 //geometry information
621 if(gelt->meshElement()==0) gelt->buildSideMeshElement();
622 MeshElement* melt = gelt->meshElement();
623 GeomMapData* mapdata=melt->geomMapData_p;
624 if(mapdata==0) mapdata=new GeomMapData(melt);
625 mapdata->computeJacobianMatrix(std::vector<real_t>(elt->dim(),0.));
626 mapdata->invertJacobianMatrix();
627 mapdata->computeDifferentialElement();
628 melt->geomMapData_p = mapdata;
629 if(nor) mapdata->computeOutwardNormal();
630
631 //shapevalues and physical points for all single quadrature
632 RefElement* relt = elt->refElt_p; //reference element related to element
633 std::set<Quadrature*>::iterator itq=quads.begin();
634 for(; itq!=quads.end(); ++itq)
635 {
636 std::map<Quadrature*,std::vector<ShapeValues> >::iterator itmq=relt->qshvs_.find(*itq);
637 dimen_t dq = (*itq)->dim(); //dimension of quadrature points
638 number_t nq = (*itq)->numberOfPoints(); //number of quadrature points
639 if(itmq==relt->qshvs_.end()) //shapevalues not computed
640 {
641 relt->qshvs_[*itq]=std::vector<ShapeValues>(nq); //to store all shapevalues at quadrature points
642 std::vector<real_t>::const_iterator itp = (*itq)->point();
643 std::vector<ShapeValues>::iterator itshv = relt->qshvs_[*itq].begin();
644 for(number_t i = 0; i < nq; i++, itp += dq, itshv++)
645 {
646 relt->computeShapeValues(itp, der);
647 *itshv = relt->shapeValues;
648 if(nbc>1) itshv->extendToVector(nbc); //extend scalar shape functions to nbc vector shape functions
649 }
650 }
651 //compute map of quadrature points
652 std::map<Quadrature*,std::vector<Point> >::iterator itmqphy= mapdata->phyPoints.find(*itq);
653 if(itmqphy==mapdata->phyPoints.end()) //create physical points
654 {
655 mapdata->phyPoints[*itq] = std::vector<Point>(nq);
656 std::vector<Point>::iterator itphy=mapdata->phyPoints[*itq].begin();
657 std::vector<real_t>::const_iterator itp = (*itq)->point();
658 for(number_t i = 0; i < nq; i++, itp += dq, ++itphy)
659 *itphy = mapdata->geomMap(itp); //quadrature point in physical space
660 }
661 }
662 }
663 }
664
665 //===================================================================================================================================
666 /*! check and set quadrature pointers of IntegrationsMethods if they are not set
667
668 if intgMethods is empty
669 if kernel is singular (Kernel::singularType!=_notSingular) add if possible
670 - a well suited singular quadrature method applied to adjacent elements (dist=0)
671 - some regular quadrature products applied to non adjacent elements (one for 0<dist<=1, one for 1<dist<=2 and one for dist>2)
672 if kernel is regular (Kernel::singularType==_notSingular) add
673 - a regular quadrature products
674 if intgMethods contains only a singular quadrature
675 if kernel is regular stop on error
676 if kernel is singular, check if singular quadrature is consistent
677 if not stop on error
678 else add some regular quadrature products applied to non adjacent elements (one for 0<dist<=1, one for 1<dist<=2 and one for dist>2)
679 if intgMethods contains only a regular quadrature
680 if kernel is regular check if it is consistent
681 if kernel is singular add a well suited singular quadrature method applied to adjacent elements (dist=0)
682 if intgMethods are more than one integration method
683 if kernel is singular,
684 check if there is only one consistent singular quadrature for dist=0, if not stop on error
685 if there is no regular quadrature, add some
686 if kernel is regular, check that all are regular quadratures
687
688 At the end of the process, produce an intgMethods vector consistent with the Kernel
689 In a second step, update pointers of quadrature if required
690 */
691 // Implemented here due to a cyclic dependance library
692 //===================================================================================================================================
setIntegrationMethods()693 void DoubleIntgBilinearForm::setIntegrationMethods()
694 {
695 const Kernel* kerp = kopus_p->kernel();
696 string_t mes="in DoubleIntgBilinearForm::setIntegrationMethods(), ";
697 const std::set<ShapeType>& shapes_u = domainu_p->meshDomain()->shapeTypes;
698 const std::set<ShapeType>& shapes_v = domainv_p->meshDomain()->shapeTypes;
699 if (shapes_u.size()!=1 || shapes_v.size()!=1)
700 {
701 where("DoubleIntgBilinearForm::setIntegrationMethods()");
702 error("domain_multiple_shapes");
703 }
704 ShapeType shu = *shapes_u.begin(), shv = *shapes_v.begin();
705
706 //analyze given integration methods
707 SingularityType st=_notsingular;
708 real_t so = 0;
709 bool hasSingPart = false, hasRegPart = false;
710 string_t sn ="Id", ssn="Id";
711 number_t nbsq=0, nbrq=0, i=0;
712 bool waitRegularPart = false;
713 if (kerp!=0) //if kerp not allocated, assume ker=Id
714 {
715 st=kerp->singularType;
716 so=kerp->singularOrder;
717 hasSingPart = kerp->singPart!=0;
718 hasRegPart = kerp->regPart!=0;
719 sn=kerp->shortname;
720 if (hasSingPart) ssn = kerp->singPart->shortname;
721 }
722 std::multimap<real_t, int_t> ims; //intg methods number indexed by distance (negative number for singular quadrature method)
723 std::multimap<real_t, int_t>::iterator itm, itme;
724 std::vector<IntgMeth>::iterator iti=intgMethods.begin();
725 for (; iti!=intgMethods.end(); ++iti, ++i) //loop on integration methods
726 {
727 const IntegrationMethod* im =iti->intgMeth;
728 if (im->imType==_quadratureIM) //single quadrature, build the product of quadrature
729 {
730 const QuadratureIM* qim = static_cast<const QuadratureIM*>(im);
731 iti->intgMeth = new ProductIM(qim->quadRule,qim->degree);
732 im = iti->intgMeth;
733 }
734 if (!im->isDoubleIM())
735 {
736 where("DoubleIntgBilinearForm::setIntegrationMethods()");
737 error("im_not_single");
738 }
739 FunctionPart pai=iti->functionPart;
740 if (pai==_regularPart && !hasRegPart)
741 {
742 where("DoubleIntgBilinearForm::setIntegrationMethods()");
743 error("null_pointer","kerp->regPart");
744 }
745 if (pai==_singularPart && !hasSingPart)
746 {
747 where("DoubleIntgBilinearForm::setIntegrationMethods()");
748 error("null_pointer","kerp->singPart");
749 }
750 SingularityType sti=im->singularType;
751 real_t soi=im->singularOrder;
752 real_t di = iti->bound;
753 if (sti ==_notsingular) //regular quadrature
754 {
755 if (st!=_notsingular && pai!=_regularPart && di==0)
756 {
757 where("DoubleIntgBilinearForm::setIntegrationMethods()");
758 error("im_not_handled",im->name);
759 }
760 ims.insert(std::pair<real_t,number_t>(di,i+1));
761 nbrq++; //regular quadrature is valid
762 }
763 else //singular quadrature
764 {
765 if (sti!=st || soi!=so) //check singularity consitency
766 {
767 where("DoubleIntgBilinearForm::setIntegrationMethods()");
768 error("im_not_handled",im->name);
769 }
770 string_t nai=im->kerName;
771 string_t na= sn;
772 if (pai==_singularPart) {na=ssn; waitRegularPart = true;}
773 if (nai!="" && nai!=na) // check if intg method is dedicated to particular kernel
774 {
775 where("DoubleIntgBilinearForm::setIntegrationMethods()");
776 error("im_not_handled",im->name);
777 }
778 if (di>0)
779 {
780 where("DoubleIntgBilinearForm::setIntegrationMethods()");
781 error("im_not_handled",im->name);
782 }
783 IntegrationMethodType imt=im->type();
784 string_t mes2="";
785 if (imt==_SauterSchwabIM && (shu!=_triangle || shv!=_triangle))
786 {
787 where("DoubleIntgBilinearForm::setIntegrationMethods()");
788 error("im_shape_only", "Sauter-Schwab", words("shape", _triangle), words("shape",shu), words("shape",shv));
789 }
790 if (imt==_DuffyIM && (shu!=_segment || shv!=_segment))
791 {
792 where("DoubleIntgBilinearForm::setIntegrationMethods()");
793 error("im_shape_only", "Duffy", words("shape", _segment), words("shape",shu), words("shape",shv));
794 }
795 if (imt==_LenoirSalles3dIM && (shu!=_triangle || shv!=_triangle))
796 {
797 where("DoubleIntgBilinearForm::setIntegrationMethods()");
798 error("im_shape_only", "Lenoir-Salles3D", words("shape", _triangle), words("shape",shu), words("shape",shv));
799 }
800 if (imt==_LenoirSalles2dIM && (shu!=_segment || shv!=_segment))
801 {
802 where("DoubleIntgBilinearForm::setIntegrationMethods()");
803 error("im_shape_only", "Lenoir-Salles2D", words("shape", _segment), words("shape",shu), words("shape",shv));
804 }
805 ims.insert(std::pair<real_t,number_t>(di,-(i+1)));
806 nbsq++; //singular quadrature is valid
807 }
808 }
809
810 // finalize analyze and add omitted integration methods by using default ones
811 if (st!=_notsingular && so<0 && domainu_p==domainv_p) //kernel is singular, nbsq should be 1
812 {
813 if (nbsq>1) //too much singular quadrature
814 {
815 where("DoubleIntgBilinearForm::setIntegrationMethods()");
816 string_t mes="";
817 itm=ims.lower_bound(0.); itme = ims.upper_bound(0.);
818 for (; itm!=itme; ++itm)
819 if (itme->second<0) mes+=" "+intgMethods[std::abs(itme->second)-1].intgMeth->name;
820 error("multiple_im",mes);
821 }
822 if (nbsq==0) //add default singular method (NOT WORK FOR MOMENT BECAUSE DuffyIMNew and SauterSchwabNew are not included, lib dependance !)
823 {
824 if (shu==_segment && shv==_segment) //use Duffy order 6
825 intgMethods.push_back(IntgMeth(DuffyIM(),_allFunction, 0.));
826 else
827 {
828 if (shu==_triangle && shv==_triangle) //use Sauter-Schwab order 3
829 intgMethods.push_back(IntgMeth(SauterSchwabIM(),_allFunction, 0.));
830 else
831 {
832 where("DoubleIntgBilinearForm::setIntegrationMethods()");
833 error("im_not_found");
834 }
835 }
836 nbsq=1;
837 }
838 }
839
840 //deal with regular quadratures
841 if (st!=_notsingular && so<0) //kernel is singular
842 {
843 if (waitRegularPart) //find if there is a regular quadrature for regularPart and dist=0
844 {
845 itm=ims.lower_bound(0.); itme = ims.upper_bound(0.);
846 for(; itm!=itme && waitRegularPart ; ++itm)
847 if(itme->second>0 && intgMethods[itme->second-1].functionPart==_regularPart) waitRegularPart=false;
848 }
849 if (waitRegularPart) // create a quadrature for regular part
850 intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,3),_regularPart, 0.));
851 if (nbrq==0) //add additional regular quadrature
852 {
853 intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,5),_allFunction, 1.));
854 intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,4),_allFunction, 2.));
855 intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,3),_allFunction, theRealMax));
856 }
857 }
858 else //kernel is regular
859 {
860 if(nbrq==0) intgMethods.push_back(IntgMeth(ProductIM(symmetrical_Gauss,3),_allFunction, theRealMax));
861 }
862
863 // update QuadratureIM pointers of regular quadratures
864 for (iti=intgMethods.begin(); iti!=intgMethods.end(); ++iti)
865 {
866 const IntegrationMethod* im =iti->intgMeth;
867 if (im->type()==_productIM) //update single quadrature pointers
868 {
869 const ProductIM* pim=static_cast<const ProductIM*>(im);
870 if (pim->getxIM()->type()==_quadratureIM) static_cast<QuadratureIM*>(pim->getxIM())->setQuadratures(shapes_v);
871 if (pim->getyIM()->type()==_quadratureIM) static_cast<QuadratureIM*>(pim->getyIM())->setQuadratures(shapes_u);
872 pim->name = pim->getxIM()->name+" x "+pim->getyIM()->name;
873 }
874 }
875 }
876
setHMIntegrationMethods()877 void DoubleIntgBilinearForm::setHMIntegrationMethods()
878 {
879 if (intgMethod_p==0) return; //no integration method
880 if (intgMethod_p->type()!=_HMatrixIM) return; //not a HMatrixIM method
881 const HMatrixIM* him= static_cast<const HMatrixIM*>(intgMethod_p);
882 intgMethods=him->intgMethods;
883 setIntegrationMethods();
884 }
885
886 // IntegrationMethods constructors
887 // implemented here due to cyclic dependance library
IntegrationMethods(IntegrationMethodType imt)888 IntegrationMethods::IntegrationMethods(IntegrationMethodType imt)
889 {
890 switch (imt)
891 {
892 case _SauterSchwabIM : add(SauterSchwabIM(),_allFunction, theRealMax); break;
893 case _DuffyIM : add(DuffyIM(),_allFunction,0.); break;
894 case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,theRealMax); break;
895 case _LenoirSalles3dIM : add(LenoirSalles3dIM(),_allFunction,theRealMax); break;
896 case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,theRealMax); break;
897 case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,theRealMax); break;
898 default:
899 where("IntegrationMethods::IntegrationMethods(IntegrationMethodType)");
900 error("im_not_handled", words("imtype", imt));
901 }
902 }
903
IntegrationMethods(QuadRule qr,number_t ord)904 IntegrationMethods::IntegrationMethods(QuadRule qr, number_t ord)
905 {
906 add(QuadratureIM(qr,ord),_allFunction,theRealMax);
907 }
908
IntegrationMethods(IntegrationMethodType imt,QuadRule qr2,number_t ord2)909 IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, QuadRule qr2, number_t ord2)
910 {
911 switch (imt)
912 {
913 case _SauterSchwabIM : add(SauterSchwabIM(),_allFunction, 0.); break;
914 case _DuffyIM : add(DuffyIM(),_allFunction,0.); break;
915 case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,0.); break;
916 case _LenoirSalles3dIM : add(LenoirSalles3dIM(),_allFunction,0.); break;
917 case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,0.); break;
918 case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,0.); break;
919 default:
920 where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, QuadRule, Number)");
921 error("im_not_handled", words("imtype", imt));
922 }
923 add(QuadratureIM(qr2,ord2),_allFunction,theRealMax);
924 }
925
IntegrationMethods(IntegrationMethodType imt,number_t ord1,real_t b1,QuadRule qr2,number_t ord2)926 IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, number_t ord1, real_t b1, QuadRule qr2, number_t ord2)
927 {
928 switch (imt)
929 {
930 case _SauterSchwabIM : add(SauterSchwabIM(ord1),_allFunction,b1); break;
931 case _DuffyIM : add(DuffyIM(ord1),_allFunction,b1); break;
932 case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
933 case _LenoirSalles3dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
934 case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,b1); break;
935 case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,b1); break;
936 default:
937 where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, Number, Real, QuadRule, Number)");
938 error("im_not_handled", words("imtype", imt));
939 }
940 add(QuadratureIM(qr2,ord2),_allFunction, theRealMax);
941 }
942
IntegrationMethods(IntegrationMethodType imt,number_t ord1,real_t b1,QuadRule qr2,number_t ord2,real_t b2,QuadRule qr3,number_t ord3)943 IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, number_t ord1, real_t b1,
944 QuadRule qr2, number_t ord2, real_t b2, QuadRule qr3, number_t ord3)
945 {
946 switch (imt)
947 {
948 case _SauterSchwabIM : add(SauterSchwabIM(ord1),_allFunction,b1); break;
949 case _DuffyIM : add(DuffyIM(ord1),_allFunction,b1); break;
950 case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
951 case _LenoirSalles3dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
952 case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,b1); break;
953 case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,b1); break;
954 default:
955 where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, Number, Real, QuadRule, Number, Real, QuadRule, Number)");
956 error("im_not_handled", words("imtype", imt));
957 }
958 add(QuadratureIM(qr2,ord2),_allFunction,b2);
959 add(QuadratureIM(qr3,ord3),_allFunction,theRealMax);
960 }
961
IntegrationMethods(IntegrationMethodType imt,number_t ord1,real_t b1,QuadRule qr2,number_t ord2,real_t b2,QuadRule qr3,number_t ord3,real_t b3,QuadRule qr4,number_t ord4)962 IntegrationMethods::IntegrationMethods(IntegrationMethodType imt, number_t ord1, real_t b1,
963 QuadRule qr2, number_t ord2, real_t b2, QuadRule qr3, number_t ord3, real_t b3, QuadRule qr4, number_t ord4)
964 {
965 switch (imt)
966 {
967 case _SauterSchwabIM : add(SauterSchwabIM(ord1),_allFunction,b1); break;
968 case _DuffyIM : add(DuffyIM(ord1),_allFunction,b1); break;
969 case _LenoirSalles2dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
970 case _LenoirSalles3dIM : add(LenoirSalles2dIM(),_allFunction,b1); break;
971 case _LenoirSalles2dIR : add(LenoirSalles2dIR(),_allFunction,b1); break;
972 case _LenoirSalles3dIR : add(LenoirSalles3dIR(),_allFunction,b1); break;
973 default:
974 where("IntegrationMethods::IntegrationMethods(IntegrationMethodType, Number, Real, QuadRule, Number, Real, QuadRule, Number, Real, QuadRule, Number)");
975 error("im_not_handled", words("imtype", imt));
976 }
977 add(QuadratureIM(qr2,ord2),_allFunction,b2);
978 add(QuadratureIM(qr3,ord3),_allFunction,b3);
979 add(QuadratureIM(qr4,ord4),_allFunction,theRealMax);
980 }
981
IntegrationMethods(QuadRule qr1,number_t ord1,real_t b1,QuadRule qr2,number_t ord2)982 IntegrationMethods::IntegrationMethods(QuadRule qr1, number_t ord1, real_t b1, QuadRule qr2, number_t ord2)
983 {
984 add(QuadratureIM(qr1,ord1),_allFunction,b1);
985 add(QuadratureIM(qr2,ord2),_allFunction,theRealMax);
986 }
987
IntegrationMethods(QuadRule qr1,number_t ord1,real_t b1,QuadRule qr2,number_t ord2,real_t b2,QuadRule qr3,number_t ord3)988 IntegrationMethods::IntegrationMethods(QuadRule qr1, number_t ord1, real_t b1, QuadRule qr2, number_t ord2, real_t b2, QuadRule qr3, number_t ord3)
989 {
990 add(QuadratureIM(qr1,ord1),_allFunction,b1);
991 add(QuadratureIM(qr2,ord2),_allFunction,b2);
992 add(QuadratureIM(qr3,ord3),_allFunction,theRealMax);
993 }
994
995 } // end of namespace xlifepp
996