1 // ********** DO NOT REMOVE THIS BANNER **********
2 // ORIG-DATE: Jan 2008
3 // -*- Mode : c++ -*-
4 //
5 // SUMMARY : Generic Fiinite Element 1d, 2d, 3d
6 // USAGE : LGPL
7 // ORG : LJLL Universite Pierre et Marie Curi, Paris, FRANCE
8 // AUTHOR : Frederic Hecht
9 // E-MAIL : frederic.hecht@ann.jussieu.fr
10 //
11
12 /*
13
14 This file is part of Freefem++
15
16 Freefem++ is free software; you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation; either version 2.1 of the License, or
19 (at your option) any later version.
20
21 Freefem++ is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with Freefem++; if not, write to the Free Software
28 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
29
30 Thank to the ARN () FF2A3 grant
31 ref:ANR-07-CIS7-002-01
32 */
33
34
35 #include <map>
36
37 #include "ufunction.hpp"
38
39 #include "error.hpp"
40 #include "RNM.hpp"
41
42 #include "Mesh3dn.hpp"
43 #include "Mesh2dn.hpp"
44
45 #include "FESpacen.hpp"
46
47 #include "splitsimplex.hpp"
48 int UniqueffId::count=0;
49 namespace Fem2D {
50
51 //template<class Element>
nbdf_d(const int ndfitem[4],const int nd[4])52 int nbdf_d(const int ndfitem[4],const int nd[4])
53 {
54 const int ndf = ndfitem[0]*nd[0] + ndfitem[1]*nd[1]+ ndfitem[2]*nd[2] + ndfitem[3]*nd[3];
55 return ndf;
56 }
57
58 //template<class Element>
nbnode_d(const int ndfitem[4],const int nd[4])59 int nbnode_d(const int ndfitem[4],const int nd[4])
60 {
61 // const int nd[]= {Element::nv, Element::ne,Element::nf,Element::nt};
62 const int ndf = nd[0]*(ndfitem[0]!=0) + nd[1]*(ndfitem[1]!=0)+ nd[2]*(ndfitem[2]!=0) + nd[3]*(ndfitem[3]!=0);
63 return ndf;
64 }
65
66 //template<class Element>
builddata_d(const int ndfitem[4],const int nd[4],int N)67 int *builddata_d(const int ndfitem[4],const int nd[4],int N)
68 {
69 // const int d=Element::Rd::d;
70 // const int nwhat=Element::nitem;
71 // const int nd[]= {Element::nv, Element::ne,Element::nf,Element::nt};
72 // const int nitem=nd[0]+nd[1]+nd[2]+nd[3];
73 // cout << " nitem="<< nitem<< endl;
74 const int ndf = nbdf_d(ndfitem,nd);
75 const int nnode=nbnode_d(ndfitem,nd);
76 int lgdata= ndf*5+N;
77 int * data = new int[lgdata];
78 int p=0;
79 for(int i=0,nw=0;i<=3;++i)
80 for(int j=0;j<nd[i];++j,nw++)// pour les what (support de node)
81 for(int k=0;k<ndfitem[i];++k)//
82 data[p++] = nw;
83 // cout << p << " " <<ndfitem[0]<< ndfitem[1]<<ndfitem[2]<<ndfitem[3]<< " " << ndf << endl;
84 assert(p==ndf);
85 for(int i=0,nw=0;i<=3;++i)
86 for(int j=0;j<nd[i];++j,nw++)// pour les what (support de node)
87 for(int k=0;k<ndfitem[i];++k)//
88 data[p++] = k;
89 // cout << p << " " << 2*ndf << " " << nitem << endl;
90 int nn=0;
91 for(int i=0;i<=3;++i)
92 {
93 int in=ndfitem[i]?1:0;
94 for(int j=0;j<nd[i];++j,nn+=in )// pour les what (support de node)
95 for(int k=0;k<ndfitem[i];++k)//
96 data[p++] = nn;
97 }
98 // cout << p << " " << 3*ndf << " " << nitem << endl;
99 for(int i=0;i<ndf*2+N;++i)
100 data[p++] = 0;
101
102 // data[p++] = 0;
103 // data[p++] = 0;
104 //cout << p << " == " << lgdata << endl;
105 assert(p== lgdata);
106 //cout << nn << " " << nnode << endl;
107 p =0;
108
109 /*
110 for(int j=0;j<4;j++)
111 {
112 for (int i=0;i<ndf;i++)
113 cout << data[p++] << " ";
114 cout << endl;
115 }
116 */
117 assert(nn==nnode);
118 return data;
119 }
120
dataTypeOfFE(const int nitemdim[4],const int dfon[4],int NN,int nbsubdivisionn,int nb_sub_femm,bool discon)121 dataTypeOfFE::dataTypeOfFE(const int nitemdim[4],const int dfon[4],int NN,int nbsubdivisionn,int nb_sub_femm,bool discon)
122 :
123 data(builddata_d(dfon,nitemdim,NN)),
124 dataalloc(data),
125 ndfonVertex(dfon[0]),
126 ndfonEdge(dfon[1]),
127 ndfonFace(dfon[2]),
128 ndfonVolume(dfon[3]),
129 NbDoF(nbdf_d(dfon,nitemdim)),
130 NbNode(nbnode_d(dfon,nitemdim)),
131 N(NN),
132 nb_sub_fem(nb_sub_femm),
133 nbsubdivision(nbsubdivisionn),
134 discontinue(discon),
135 DFOnWhat(data+0*NbDoF),
136 DFOfNode(data+1*NbDoF),
137 NodeOfDF(data+2*NbDoF),
138 fromFE(data+3*NbDoF),
139 fromDF(data+4*NbDoF),
140 fromASubFE(data+3*NbDoF),
141 fromASubDF(data+4*NbDoF) ,
142 dim_which_sub_fem(data+5*NbDoF)
143 {}
144
145
builddata_d(const int nitemdim[4],const KN<dataTypeOfFE const * > & teb)146 int *builddata_d(const int nitemdim[4],const KN< dataTypeOfFE const *> &teb)
147 {
148 const int k = teb.N();
149 KN<int> NN(k+1), DF(k+1) , comp(k+1);
150 map< dataTypeOfFE const *,int> m;
151 int i=k,j;
152 while(i--) // on va a l'envert pour avoir comp[i] <=i
153 m[teb[i]]=i;
154 // l'ordre comp est important comp est croissant mais pas de pb.
155 i=k;
156 while(i--)
157 comp[i]=m[teb[i]]; // comp[i] <=i
158 int n=0,N=0;
159 for ( j=0;j<k;j++)
160 {NN[j]=N;N+=teb[j]->N;}
161 NN[k] = N;
162 // reservation des interval en df
163 n=0;
164 for ( j=0;j<k;j++)
165 { DF[j]=n;n+=teb[j]->NbDoF;}
166 DF[k] = n;
167
168 int NbDoF=0;
169 int dfon[4]={0,0,0,0};
170 int nbsubdivision=0;
171 int discon=0;
172 for (int i=0;i<k;++i)
173 {
174 NbDoF += teb[i]->NbDoF;
175 dfon[0] += teb[i]->ndfonVertex;
176 dfon[1] += teb[i]->ndfonEdge;
177 dfon[2] += teb[i]->ndfonFace;
178 dfon[3] += teb[i]->ndfonVolume;
179 nbsubdivision = max(nbsubdivision,teb[i]->nbsubdivision);
180 discon = discon || teb[i]->discontinue; // bof bof 1 FE discontinue => discontinue
181 }
182 int nwhat=15; // 15 = 4+6+1+1 (nb of support item (what) : vertex, edges, fqces, tet)
183
184 int ostart=nwhat;
185 int * data0=new int[ostart+7*NbDoF+N];
186 int * data=data0+ostart;
187 int * data1=data+5*NbDoF;
188
189 int c=0;
190 KN<int> w(nwhat),nn(nwhat);
191
192 w=0;
193 nn=0;
194
195
196 for ( j=0;j<k;j++)
197 for ( i=0;i<teb[j]->NbDoF;i++)
198 nn[teb[j]->DFOnWhat[i]]++;
199 int nbn=0;
200 for( j=0;j<nwhat;j++)
201 if (nn[j]) nn[j]=nbn++;
202 else nn[j]=-1;
203 KN<int> dln(nwhat);
204 dln=0;
205 // nn donne numero de noeud sur what
206 for ( j=0;j<k;j++)
207 for ( i=0;i<teb[j]->NbDoF;i++)
208 data[c++] = teb[j]->DFOnWhat[i];
209
210 for ( j=0;j<k;j++)
211 {
212 int cc=c;
213 for ( i=0;i<teb[j]->NbDoF;i++)
214 data[c++] = teb[j]->DFOfNode[i]+dln[teb[j]->DFOnWhat[i]];
215 for ( i=0;i<teb[j]->NbDoF;i++)
216 dln[teb[j]->DFOnWhat[i]]=Max(dln[teb[j]->DFOnWhat[i]],data[cc++]+1);
217 }
218
219
220 for ( j=0;j<k;j++)
221 {
222 // w renumerotation des noeuds
223 // Ok si un noeud par what
224 for ( i=0;i<teb[j]->NbDoF;i++)
225 data[c++] = nn[teb[j]->DFOnWhat[i]];
226 }
227
228 for ( j=0;j<k;j++)
229 for ( i=0;i<teb[j]->NbDoF;i++)
230 data[c++] = j; // node from of FE
231
232
233 for ( j=0;j<k;j++)
234 for ( i=0;i<teb[j]->NbDoF;i++)
235 data[c++] = i; // node from of df in FE
236 // error -- here
237 //in case of [P2,P2],P1
238 // we expect 0,0,1 and we get 0 1 2
239 // => wrong BC ????
240 c+=2*n; // on saute le deux tableau en plus (cf data1.)
241
242
243 int xx=0;
244 for (j=0;j<k;j++)
245 {
246 int xxx=xx;
247 for (i=0;i<teb[j]->N;i++)
248 {
249 data[c] = teb[j]->dim_which_sub_fem[i]+xx;
250 xxx=Max(xxx,data[c]+1);
251 c++;
252 }
253 xx=xxx;
254 }
255
256
257 // ou dans la partie miminal element finite atomic
258
259 int ci=n;
260 int cj=0;
261 int ccc=0;
262 for ( j=0;j<k;ccc+=teb[j++]->nb_sub_fem)
263
264 for ( i=0;i<teb[j]->NbDoF;i++)
265 {
266 int il= teb[j]->fromASubDF[i];
267 int jl= teb[j]->fromASubFE[i];
268 data1[ci++]=il;
269 data1[cj++]=ccc+jl;
270 }
271
272 int nb_sub_fem=ccc;
273
274 ffassert(c== 7*n+N);
275 /* int cc=0;
276 cout << " Data : " << endl;
277 for ( i=0;i<5;i++) {
278 for (j=0;j<n;j++)
279 cout << " " << data[cc++];
280 cout << endl;}
281 cout << " which " ;
282 for (i=0;i<N;i++)
283 cout << " " << data[cc++];
284 cout << endl;*/
285
286
287 for(int i=0;i<4;++i)
288 data0[i]=dfon[i];
289
290 data0[4]=NbDoF;
291 data0[5]=nbn;// NbNode
292 data0[6]=N;
293 data0[7]=nb_sub_fem;
294 data0[8]=nbsubdivision;
295 data0[9]=discon;
296
297 return data0;
298 }
299
300
dataTypeOfFE(const int nitemdim[4],const KN<dataTypeOfFE const * > & tef)301 dataTypeOfFE::dataTypeOfFE(const int nitemdim[4],const KN< dataTypeOfFE const *> & tef)
302 :
303 data(builddata_d(nitemdim,tef)),
304 dataalloc(data),
305 ndfonVertex(data[0]),
306 ndfonEdge(data[1]),
307 ndfonFace(data[2]),
308 ndfonVolume(data[3]),
309 NbDoF(data[4]),
310 NbNode(data[5]),
311 N(data[6]),
312 nb_sub_fem(data[7]),
313 nbsubdivision(data[8]),
314 discontinue(data[9]),
315 DFOnWhat(data+15+0*NbDoF),
316 DFOfNode(data+15+1*NbDoF),
317 NodeOfDF(data+15+2*NbDoF),
318 fromFE(data+15+3*NbDoF),
319 fromDF(data+15+4*NbDoF),
320 fromASubFE(data+15+5*NbDoF),
321 fromASubDF(data+15+6*NbDoF) ,
322 dim_which_sub_fem(data+15+7*NbDoF)
323 {}
324
325 template<class Mesh>
init(InterpolationMatrix<RdHat> & M,FElement * pK,int odf,int ocomp,int * pp) const326 void GTypeOfFESum<Mesh>::init(InterpolationMatrix<RdHat> & M,FElement * pK,int odf,int ocomp,int *pp) const
327 {
328 // a faire ..... cas matrix invariante
329 assert(0);
330 }
331
332 template<class Mesh>
GTypeOfFESum(const KN<GTypeOfFE<Mesh> const * > & t)333 GTypeOfFESum<Mesh>::GTypeOfFESum(const KN< GTypeOfFE<Mesh> const *> & t)
334 :
335 GTypeOfFE<Mesh>(t),
336 k(t.N()),
337 teb(t),
338 NN(k+1),
339 DF(k+1) ,
340 comp(k+1) {Build();}
341
342 template<class Mesh>
kn(const GFESpace<Mesh> ** tt,int kk)343 static KN< GTypeOfFE<Mesh> const *> kn(const GFESpace<Mesh> ** tt,int kk)
344 {
345 KN< GTypeOfFE<Mesh> const *> r(kk);
346 for(int i=0;i<kk;++i)
347 { r[i]=tt[i]->TFE[0];ffassert(tt[i]->TFE.constant());}
348 return r;
349 }
350 template<class Mesh>
kn(const GFESpace<Mesh> & tt,int kk)351 static KN< GTypeOfFE<Mesh> const *> kn(const GFESpace<Mesh> & tt,int kk)
352 {
353 return KN< GTypeOfFE<Mesh> const *> (kk,tt.TFE[0]);
354 }
355
356 template<class Mesh>
GTypeOfFESum(const GFESpace<Mesh> ** tt,int kk)357 GTypeOfFESum<Mesh>::GTypeOfFESum(const GFESpace<Mesh> ** tt,int kk)
358 :
359 GTypeOfFE<Mesh>(kn(tt,kk)),
360 k(kk),
361 teb(kn(tt,kk)),
362 NN(k+1),
363 DF(k+1) ,
364 comp(k+1) {Build();}
365
366 template<class Mesh>
GTypeOfFESum(const GFESpace<Mesh> & tt,int kk)367 GTypeOfFESum<Mesh>::GTypeOfFESum(const GFESpace<Mesh> & tt,int kk)
368 :
369 GTypeOfFE<Mesh>(kn(tt,kk)),
370 k(kk),
371 teb(kn(tt,kk)),
372 NN(k+1),
373 DF(k+1) ,
374 comp(k+1) {Build();}
375
376 template<class Mesh>
Build()377 void GTypeOfFESum<Mesh>::Build()
378 {
379 bool debug=verbosity>5;;
380 {
381 const KN< GTypeOfFE<Mesh> const *> & t=teb;
382 map<const GTypeOfFE<Mesh> *,int> m;
383 int i=k,j;
384 while(i--) // on va a l'envert pour avoir comp[i] <=i
385 m[teb[i]]=i;
386 // l'ordre comp est important comp est croissant mais pas de pb.
387 i=k;
388 while(i--)
389 comp[i]=m[teb[i]]; // comp[i] <=i
390
391 // reservatition des intervalles en espaces
392 int n=0,N=0;
393 for ( j=0;j<k;j++)
394 {NN[j]=N;N+=teb[j]->N;}
395 NN[k] = N;
396 // reservation des interval en df
397 n=0;
398 for ( j=0;j<k;j++)
399 { DF[j]=n;n+=teb[j]->NbDoF;}
400 DF[k] = n;
401 }
402 int ii=0;
403 for (int i=0;i<k;++i)
404 {
405 for (int j=0;j<teb[i]->nb_sub_fem;++j)
406 this->Sub_ToFE[ii++]=teb[i]->Sub_ToFE[j];
407 }
408 assert(ii==this->nb_sub_fem );
409
410 int c=0,c0=0, fcom=0;
411 for (int i=0;i<this->nb_sub_fem;i++)
412 {
413 int N=this->Sub_ToFE[i]->N;
414 int ndofi=this->Sub_ToFE[i]->NbDoF;
415 this->first_comp[i]= fcom;
416 this->last_comp[i]= fcom+N;
417 fcom += N;
418
419 for(int j=0;j<N;++j)
420 {
421 this->begin_dfcomp[c] = c0 + this->Sub_ToFE[i]->begin_dfcomp[j] ;
422 this->end_dfcomp[c] = c0 + this->Sub_ToFE[i]->end_dfcomp[j] ;
423 c++;
424 }
425 c0+=ndofi;
426
427 }
428 if(debug)
429 {
430 cout <<" NbDoF : " << this->NbDoF <<endl;
431 for(int i=0;i<this->N;++i)
432 cout << " comp " << i << " ["<<this->begin_dfcomp[i]<<", "<< this->end_dfcomp[i]<< "[\n";
433 }
434
435 // construction de l'interpolation .
436
437 int npi=0;
438 int nci=0;
439 bool var=true;
440 for (int i=0;i<this->nb_sub_fem;i++)
441 {
442 npi +=this->Sub_ToFE[i]->NbPtforInterpolation;
443 nci +=this->Sub_ToFE[i]->NbcoefforInterpolation;
444 var = var && this->Sub_ToFE[i]->invariantinterpolationMatrix;
445 }
446 assert(this->NbcoefforInterpolation== nci);
447 this->invariantinterpolationMatrix=var;
448 // this->pInterpolation.init(nci);
449 // this->cInterpolation.init(nci);
450 // this->dofInterpolation.iniy(nci);
451 KN<int> opi(this->nb_sub_fem);// offset numumber point intgartion
452 {
453 map<RdHat,int,lessRd> mpt;
454 numPtInterpolation.init(npi);
455 int npp=0,kkk=0;
456 KN<RdHat> Ptt(npi);
457 for (int i=0;i<this->nb_sub_fem;i++)
458 {
459 opi[i]=kkk;
460 const GTypeOfFE<Mesh> &ti=*this->Sub_ToFE[i];
461
462 for(int p=0;p<ti.NbPtforInterpolation;++p,++kkk)
463 {
464 Ptt[kkk]=ti.PtInterpolation[p];
465 if( mpt.find(Ptt[kkk]) == mpt.end())
466 mpt[Ptt[kkk]]=npp++;
467 numPtInterpolation[kkk]=mpt[Ptt[kkk]];
468 if(verbosity>100)
469 cout << " p= "<< p << " [ " << Ptt[kkk]<< "] , "<< kkk<< " "<< npp<< " " << numPtInterpolation[kkk]<< endl;;
470
471 }
472 }
473 assert(this->NbPtforInterpolation==0);
474 if(verbosity>5)
475 cout << npp;
476 this->NbPtforInterpolation=npp;
477 this->PtInterpolation.init(npp);
478 for(int i=0;i<npp;++i)
479 this->PtInterpolation[numPtInterpolation[i]]=Ptt[i];
480 }
481
482 int oc=0,odof=0;
483 for (int i=0,k=0;i<this->nb_sub_fem;i++)
484 {
485 const GTypeOfFE<Mesh> &ti=*this->Sub_ToFE[i];
486 for(int j=0;j<ti.NbcoefforInterpolation; ++j,++k)
487 {
488 this->pInterpolation[k] = numPtInterpolation[opi[i]+ti.pInterpolation[j]];
489 this->cInterpolation[k] = ti.cInterpolation[j]+oc;
490 this->dofInterpolation[k] = ti.dofInterpolation[j]+odof;
491 this->coefInterpolation[k]=ti.coefInterpolation[j];
492 }
493 oc += ti.N;
494 odof += ti.NbDoF;
495 }
496 if(verbosity>100)
497 cout << " **GTypeOfFESum<Mesh>::Build() " <<this->pInterpolation <<endl;
498 assert(c==this->N);
499 }
500
501
set(const Mesh & Th,const Element & K,InterpolationMatrix<RdHat> & M,int oocoef,int oodf,int * nnump) const502 template<class Mesh> void GTypeOfFESum<Mesh>::set(const Mesh & Th,const Element & K,InterpolationMatrix<RdHat> & M,int oocoef,int oodf,int *nnump ) const
503 {
504 int op=0,oc=0,odof=oodf,ocoef=oocoef;
505 assert(nnump==0);
506 for (int i=0,k=0;i<this->nb_sub_fem;i++)
507 {
508 const GTypeOfFE<Mesh> &ti=*this->Sub_ToFE[i];
509 if(!ti.invariantinterpolationMatrix)
510 ti.set(Th,K,M,ocoef,odof,&numPtInterpolation[op]);
511 oc += ti.N;
512 odof += ti.NbDoF;
513 ocoef += ti.NbcoefforInterpolation;
514 op += ti.NbPtforInterpolation;
515
516 }
517 if( verbosity > 100) cout << " GTypeOfFESum set "<< this->coefInterpolation << endl;
518 }
519
520 template<class MMesh>
GFESpace(const GFESpace & Vh,int kk,int nbequibe,int * equibe)521 GFESpace<MMesh>::GFESpace(const GFESpace & Vh,int kk,int nbequibe,int *equibe)
522 :
523 GFESpacePtrTFE<MMesh>(new GTypeOfFESum<MMesh>(Vh,kk)),
524 DataFENodeDF(Vh.Th.BuildDFNumbering(this->ptrTFE->ndfonVertex,this->ptrTFE->ndfonEdge,this->ptrTFE->ndfonFace,this->ptrTFE->ndfonVolume,nbequibe,equibe)),
525 Th(Vh.Th),
526 TFE(1,0,this->ptrTFE),
527 cmesh(Th),
528 N(TFE[0]->N),
529 Nproduit(kk),
530 nb_sub_fem(TFE[0]->nb_sub_fem),
531 dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
532 maxNbPtforInterpolation(TFE[0]->NbPtforInterpolation),
533 maxNbcoefforInterpolation(TFE[0]->NbcoefforInterpolation)
534
535 {
536 }
537
538 template<class MMesh>
GFESpace(const GFESpace ** pVh,int kk,int nbequibe,int * equibe)539 GFESpace<MMesh>::GFESpace(const GFESpace ** pVh,int kk,int nbequibe,int *equibe)
540 :
541 GFESpacePtrTFE<MMesh>(new GTypeOfFESum<MMesh>(pVh,kk)),
542 DataFENodeDF((**pVh).Th.BuildDFNumbering(this->ptrTFE->ndfonVertex,this->ptrTFE->ndfonEdge,this->ptrTFE->ndfonFace,this->ptrTFE->ndfonVolume,nbequibe,equibe)),
543 Th((**pVh).Th),
544 TFE(1,0,this->ptrTFE),
545 cmesh(Th),
546 N(TFE[0]->N),
547 Nproduit(FirstDfOfNodeData ? 1 :MaxNbDFPerNode),
548 nb_sub_fem(TFE[0]->nb_sub_fem),
549 dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
550 maxNbPtforInterpolation(TFE[0]->NbPtforInterpolation),
551 maxNbcoefforInterpolation(TFE[0]->NbcoefforInterpolation)
552
553 {
554 long snbdf=0;
555 for(int i=0;i<kk;++i)
556 snbdf += pVh[i]->NbOfDF;
557 if( snbdf !=NbOfDF)
558 cerr << " Problem build of GFESpace (3d) (may be : due to periodic Boundary condition missing ) FH " << endl
559 << " The number of DF must be " << snbdf << " and it is " << NbOfDF <<endl;
560 ffassert(snbdf == NbOfDF );
561
562
563 for(int i=0;i<kk;++i)
564 ffassert(&Th==&pVh[i]->Th);
565 }
566
567 template<class MMesh>
568 template<class R>
newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<typename MMesh::RdHat> & Psub,KN<int> & Ksub,int op_U) const569 KN<R> GFESpace<MMesh>::newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<typename MMesh::RdHat> &Psub,KN<int> &Ksub,int op_U) const
570 {
571 const int d = Rd::d;
572 Rd *Ps=0;
573 int *Ks=0;
574 int nsb = TFE[0]->nbsubdivision;
575 int nvsub,nksub;
576 SplitSimplex<Rd>(nsb, nvsub, Ps, nksub , Ks);
577 ffassert( Psub.unset());
578 ffassert( Ksub.unset());
579 Psub.set(Ps,nvsub);
580 Ksub.set(Ks,nksub*(d+1));
581 lg= nvsub*Th.nt;
582 KN<R> v(lg);
583 for (int k=0,i=0;k<Th.nt;k++)
584 {
585 FElement K=(*this)[k];
586 for(int l=0;l<nvsub;l++)
587 v[i++] = K(Psub[l], U, componante, op_U) ;
588
589 }
590 return KN<R>(true,v);// to remove the copy.
591 }
592
593
594 template< >
595 template<class R>
newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<typename MeshS::RdHat> & Psub,KN<int> & Ksub,int op_U) const596 KN<R> GFESpace<MeshS>::newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<typename MeshS::RdHat> &Psub,KN<int> &Ksub,int op_U) const
597 {
598 typedef typename MeshS::RdHat RdHat;
599 const int dHat = RdHat::d;
600 MeshS::RdHat *Ps=0;
601 int *Ks=0;
602 int nsb = TFE[0]->nbsubdivision;
603 int nvsub,nksub;
604 SplitSimplex<RdHat>(nsb, nvsub, Ps, nksub , Ks);
605 ffassert( Psub.unset());
606 ffassert( Ksub.unset());
607 Psub.set(Ps,nvsub);
608 Ksub.set(Ks,nksub*(dHat+1));
609 lg= nvsub*Th.nt;
610 KN<R> v(lg);
611 for (int k=0,i=0;k<Th.nt;k++)
612 {
613 FElementS K=(*this)[k];
614 for(int l=0;l<nvsub;l++)
615 v[i++] = K(Psub[l], U, componante, op_U) ;
616
617 }
618 return KN<R>(true,v);// to remove the copy.
619 }
620
621
622 template< >
623 template<class R>
newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<typename MeshL::RdHat> & Psub,KN<int> & Ksub,int op_U) const624 KN<R> GFESpace<MeshL>::newSaveDraw(const KN_<R> & U,int componante,int & lg,KN<typename MeshL::RdHat> &Psub,KN<int> &Ksub,int op_U) const
625 {
626 typedef typename MeshL::RdHat RdHat;
627 const int dHat = RdHat::d;
628 MeshL::RdHat *Ps=0;
629 int *Ks=0;
630 int nsb = TFE[0]->nbsubdivision;
631 int nvsub,nksub;
632 SplitSimplex<RdHat>(nsb, nvsub, Ps, nksub , Ks);
633 ffassert( Psub.unset());
634 ffassert( Ksub.unset());
635 Psub.set(Ps,nvsub);
636 Ksub.set(Ks,nksub*(dHat+1));
637 lg= nvsub*Th.nt;
638 KN<R> v(lg);
639 for (int k=0,i=0;k<Th.nt;k++)
640 {
641 FElementL K=(*this)[k];
642 for(int l=0;l<nvsub;l++)
643 v[i++] = K(Psub[l], U, componante, op_U) ;
644
645 }
646 return KN<R>(true,v);// to remove the copy.
647 }
648
649
650
651 /*
652 template<class MMesh>
653 KN<double> GFESpace<MMesh>::newSaveDraw(const KN_<R> & U,int composante,int & lg,int & nsb) const
654 {
655 nsb = TFE[0]->nbsubdivision;
656 int nsbv = NbOfSubInternalVertices(nsb,d);
657 lg = nsbv*Th.nt;
658 cout << "newSaveDraw what: nt " << Th.nt << " " << nsbv << " " << lg << endl;
659 KN<double> v(lg);
660 ffassert(v);
661 for (int k=0,i=0;k<Th.nt;k++)
662 {
663 (*this)[k].SaveDraw( U,composante,&v[i]);
664 i+=nsbv;
665 }
666 return KN<double>(true,v);// to remove the copy.
667 }
668 */
669 // explicite instance..
670 template class GTypeOfFESum<Mesh2>;
671 template class GTypeOfFESum<Mesh3>;
672 template class GTypeOfFESum<MeshS>;
673 template class GTypeOfFESum<MeshL>;
674 template class GFESpace<Mesh1>;
675 template class GFESpace<Mesh2>;
676 template class GFESpace<Mesh3>;
677 template class GFESpace<MeshS>;
678 template class GFESpace<MeshL>;
679
680 template KN<double> GFESpace<MeshL>::newSaveDraw<double>(const KN_<double> & U,int componante,int & lg,KN<typename MeshL::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
681 template KN<double> GFESpace<MeshS>::newSaveDraw<double>(const KN_<double> & U,int componante,int & lg,KN<typename MeshS::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
682 template KN<double> GFESpace<Mesh3>::newSaveDraw<double>(const KN_<double> & U,int componante,int & lg,KN<typename Mesh3::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
683 template KN<double> GFESpace<Mesh2>::newSaveDraw<double>(const KN_<double> & U,int componante,int & lg,KN<typename Mesh2::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
684 template KN<double> GFESpace<Mesh1>::newSaveDraw<double>(const KN_<double> & U,int componante,int & lg,KN<typename Mesh1::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
685
686 typedef std::complex<double> Complex;
687 template KN<Complex> GFESpace<MeshL>::newSaveDraw<Complex>(const KN_<Complex> & U,int componante,int & lg,KN<typename MeshL::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
688 template KN<Complex> GFESpace<MeshS>::newSaveDraw<Complex>(const KN_<Complex> & U,int componante,int & lg,KN<typename MeshS::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
689 template KN<Complex> GFESpace<Mesh3>::newSaveDraw<Complex>(const KN_<Complex> & U,int componante,int & lg,KN<typename Mesh3::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
690 template KN<Complex> GFESpace<Mesh2>::newSaveDraw<Complex>(const KN_<Complex> & U,int componante,int & lg,KN<typename Mesh2::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
691 template KN<Complex> GFESpace<Mesh1>::newSaveDraw<Complex>(const KN_<Complex> & U,int componante,int & lg,KN<typename Mesh1::RdHat> &Psub,KN<int> &Ksub,int op_U) const ;
692
693
694 }
695