NeoPZ
PrismExtend.cpp
Go to the documentation of this file.
1 
5 #include "PrismExtend.h"
6 #include "tpzpoint.h"
7 
8 using namespace std;
9 
10 #include "pzlog.h"
11 #ifdef LOG4CXX
12 static LoggerPtr logger(Logger::getLogger("pz.topology.prismextend"));
13 #endif
14 
15 namespace pztopology {
16 
17  template<class TFather>
18  Pr<TFather>::Pr()
19  {
20  }
21 
22  template<class TFather>
23  Pr<TFather>::~Pr()
24  {
25  }
26 
27  template<class TFather>
28  int Pr<TFather>::NSideNodes(int side)
29  {
30  int ftns = side/TFather::NSides;
31  if(ftns < 2) return TFather::NSideNodes(side % TFather::NSides);
32  return TFather::NSideNodes(side % TFather::NSides) *2;
33  }
34 
35  template<class TFather>
36  int Pr<TFather>::SideNodeLocId(int side, int node)
37  {
38  int ftns = side/TFather::NSides;
39  int fatherside = side%TFather::NSides;
40  if(ftns == 0) return TFather::SideNodeLocId(side,node);
41  if(ftns == 1) return TFather::SideNodeLocId(side-TFather::NSides,node)+TFather::NCornerNodes;
42  int fatherlevel = node/TFather::NSideNodes(fatherside);
43  if(fatherlevel == 0) return TFather::SideNodeLocId(fatherside,node);
44  if(fatherlevel == 1) return TFather::SideNodeLocId(fatherside,node-TFather::NSideNodes(fatherside))+TFather::NCornerNodes;//TFather::NSideNodes(fatherside);
45  PZError << "Pr<TFather>::SideNodeLocId inconsistent side or node " << side
46  << ' ' << node << endl;
47  return -1;
48  }
49 
50  template<class TFather>
51  int Pr<TFather>::SideDimension(int side) {
52  if(side<0 || side >= NSides) {
53  PZError << "Pr<TFather>::SideDimension side " << side << endl;
54  return -1;
55  }
56  int ftns = side/TFather::NSides;
57  int fatherside = side%TFather::NSides;
58  if(ftns < 2) return TFather::SideDimension(fatherside);
59  return TFather::SideDimension(fatherside)+1;
60  }
61 
62  template<class TFather>
63  void Pr<TFather>::HigherDimensionSides(int side, TPZStack<int> &high)
64  {
65  if(side <0 || side >= NSides) {
66  PZError << "Pr<TFather>::HigherDimensionSides side "<< side << endl;
67  }
68  int ftns = side/TFather::NSides;
69  int fatherside = side%TFather::NSides;
70  int nhbefore = high.NElements();
71  TFather::HigherDimensionSides(fatherside,high);
72  int nhafter = high.NElements();
73  int is;
74  if(ftns < 2)
75  {
76  high.Push(fatherside+2*TFather::NSides);
77  for(is=nhbefore; is<nhafter; is++)
78  {
79  high.Push(high[is]+2*TFather::NSides);
80  high[is] += ftns*TFather::NSides;
81  }
82  }
83  if(ftns == 2)
84  {
85  for(is=nhbefore; is<nhafter; is++)
86  {
87  high[is] += ftns*TFather::NSides;
88  }
89  }
90  }
91 
92  template<class TFather>
93  void Pr<TFather>::CenterPoint(int side, TPZVec<REAL> &center) {
94  if (center.size()!=Dimension) {
95  DebugStop();
96  }
97  int ftns = side/TFather::NSides;
98  int fatherside = side%TFather::NSides;
99  TFather::CenterPoint(fatherside,center);
100  switch(ftns)
101  {
102  case 0:
103  center[Dimension-1] = -1.;
104  break;
105  case 1:
106  center[Dimension-1] = 1.;
107  break;
108  case 2:
109  default:
110  center[Dimension-1] = 0.;
111  }
112  }
113 
114  template<class TFather>
115  TPZTransform<> Pr<TFather>::TransformElementToSide(int side){
116 
117  if(side<0 || side>=NSides)
118  {
119  PZError << "Pr<TFather>::TransformElementToSide called with side error\n";
120  return TPZTransform<>(0,0);
121  }
122 
123  int sidedim = SideDimension(side);
124  TPZTransform<> t(sidedim,Dimension);//t(dimto,2)
125  if(side == NSides) return t;
126  int ftns = side/TFather::NSides;
127  int fatherside = side%TFather::NSides;
128  TPZTransform<> fathert = TFather::TransformElementToSide(fatherside);
129  int fathersidedimension = sidedim;
130  if(ftns == 2) fathersidedimension--;
131  int id,jd;
132  for(id=0; id<fathersidedimension; id++)
133  {
134  t.Sum()(id,0) = fathert.Sum()(id,0);
135  for(jd=0; jd<TFather::Dimension; jd++)
136  {
137  t.Mult()(id,jd) = fathert.Mult()(id,jd);
138  }
139  }
140  if(ftns == 2)
141  {
142  t.Mult()(sidedim-1,Dimension-1) = 1.;
143  }
144  return t;
145  }
146 
147  template<class TFather>
148  TPZTransform<> Pr<TFather>::TransformSideToElement(int side){
149 
150  if(side<0 || side>=NSides){
151  PZError << "Pr<TFather>::TransformSideToElement side out range\n";
152  return TPZTransform<>(0,0);
153  }
154  int sidedim = SideDimension(side);
155  TPZTransform<> t(Dimension,sidedim);
156  if(Dimension == sidedim) return t;
157  int ftns = side/TFather::NSides;
158  int fatherside = side%TFather::NSides;
159  TPZTransform<> fathert = TFather::TransformSideToElement(fatherside);
160  int fathersidedimension = sidedim;
161  if(ftns == 2) fathersidedimension--;
162  int id,jd;
163  for(id=0; id<TFather::Dimension; id++)
164  {
165  t.Sum()(id,0) = fathert.Sum()(id,0);
166  for(jd=0; jd<fathersidedimension; jd++)
167  {
168  t.Mult()(id,jd) = fathert.Mult()(id,jd);
169  }
170  }
171  switch(ftns)
172  {
173  case 0:
174  t.Sum()(Dimension-1,0) = -1.;
175  break;
176  case 1:
177  t.Sum()(Dimension-1,0) = 1.;
178  break;
179  case 2:
180  t.Mult()(Dimension-1,sidedim-1) = 1.;
181  break;
182  }
183  return t;
184 
185  }
186 
187  template<class TFather>
188  TPZTransform<> Pr<TFather>::SideToSideTransform(int sidefrom, int sideto)
189  {
190  if(sidefrom <0 || sidefrom >= NSides || sideto <0 || sideto >= NSides) {
191  PZError << "Pr<TFather>::HigherDimensionSides sidefrom "<< sidefrom <<
192  ' ' << sideto << endl;
193  return TPZTransform<>(0);
194  }
195  if(sidefrom == sideto) {
196  return TPZTransform<>(SideDimension(sidefrom));
197  }
198  if(sidefrom == NSides-1) {
199  return TransformElementToSide(sideto);
200  }
201  TPZStack<int> highsides;
202  HigherDimensionSides(sidefrom,highsides);
203  int nhigh = highsides.NElements();
204  int is;
205  for(is=0; is<nhigh; is++) {
206  if(highsides[is] == sideto) {
207  int dfr = SideDimension(sidefrom);
208  int dto = SideDimension(sideto);
209  TPZTransform<> trans(dto,dfr),toelement(Dimension,sidefrom),toside(sideto,Dimension);
210  toelement = TransformSideToElement(sidefrom);
211  toside = TransformElementToSide(sideto);
212  trans = toside.Multiply(toelement);
213  return trans;
214  }
215  }
216  PZError << "Pr<TFather>::SideToSideTransform highside not found sidefrom "
217  << sidefrom << ' ' << sideto << endl;
218  return TPZTransform<>(0);
219  }
220 
221  template<class TFather>
222  int Pr<TFather>::NumSides() {
223  return TFather::NumSides()*3;
224  }
225 
226 
227  template<class TFather>
228  int Pr<TFather>::NContainedSides(int side) {
229  int ftns = side/TFather::NSides;
230  int fatherside = side%TFather::NSides;
231  if(ftns<2) return TFather::NContainedSides(fatherside);
232  return TFather::NContainedSides(fatherside)*3;
233  }
234 
235  template<class TFather>
236  int Pr<TFather>::ContainedSideLocId(int side,int c) {
237  int nsconnect = NContainedSides(side);
238  if(c >= nsconnect)
239  {
240  PZError << "Pr<TFather>::ContainedSideLocId, side = " << side << " connect = " << c << endl;
241  return -1;
242  }
243  int ftns = side/TFather::NSides;
244  int fatherside = side%TFather::NSides;
245  if(ftns == 0) return TFather::ContainedSideLocId(fatherside,c);
246  if(ftns == 1) return TFather::ContainedSideLocId(fatherside,c)+TFather::NSides;
247  int nsideconfather = TFather::NContainedSides(fatherside);
248  int confather = c%nsideconfather;
249  int level = c/nsideconfather;
250  return TFather::ContainedSideLocId(fatherside,confather)+level*TFather::NSides;
251  }
252 
253  template<class TFather>
254  void Pr<TFather>::LowerDimensionSides(int side,TPZStack<int> &smallsides)
255  {
256  smallsides.Resize(0);
257  int nsidecon = NContainedSides(side);
258  int is;
259  for(is=0; is<nsidecon-1; is++)
260  smallsides.Push(ContainedSideLocId(side,is));
261  }
262 
263  template<class TFather>
264  void Pr<TFather>::LowerDimensionSides(int side,TPZStack<int> &smallsides, int DimTarget)
265  {
266  smallsides.Resize(0);
267  int nsidecon = NContainedSides(side);
268  for(int is = 0; is < nsidecon - 1; is++) {
269  if (SideDimension(ContainedSideLocId(side,is)) == DimTarget) smallsides.Push(ContainedSideLocId(side,is));
270  }
271  }
272 
274  template<class TFather>
275  REAL Pr<TFather>::RefElVolume()
276  {
277  return 2.0*TFather::RefElVolume();
278  }
279 
280  template<class TFather>
281  TPZIntPoints *Pr<TFather>::CreateSideIntegrationRule(int side, int order) {
282 
283  TPZIntPoints *integ = TFather::CreateSideIntegrationRule(side%TFather::NSides,order);
284  if(side< 2*TFather::NSides) {
285  return integ;
286  }
287  TPZIntPoints *integext = integ->PrismExtend(order);
288  delete integ;
289  return integext;
290 
291  }
292 
293  template<class TFather>
294  std::string Pr<TFather>::StrType()
295  {
296  return "Pr:" + MElementType_Name(TFather::Type());
297  }
298 
299  template<class TFather>
300  std::string Pr<TFather>::StrType(int side)
301  {
302  if(side < 2*TFather::NSides)
303  {
304  return MElementType_Name(TFather::Type(side%(TFather::NSides)));
305  }
306  return "Pr:" + MElementType_Name(TFather::Type(side%(TFather::NSides)));
307  }
308 
309  template<class TFather>
310  void Pr<TFather>::MapToSide(int side, TPZVec<REAL> &InternalPar, TPZVec<REAL> &SidePar, TPZFMatrix<REAL> &JacToSide) {
311  TPZTransform<> Transf = SideToSideTransform(NSides - 1, side);
312  SidePar.Resize(SideDimension(side));
313  Transf.Apply(InternalPar,SidePar);
314 
315  JacToSide = Transf.Mult();
316  }
317 
321  template<class TFather>
322  void Pr<TFather>::Diagnostic()
323  {
324 #ifdef LOG4CXX
325  if (logger->isDebugEnabled())
326  {
327  std::stringstream sout;
328  sout << __PRETTY_FUNCTION__ ;
329  sout << " Dimension : " << Dimension << " NCornerNodes " << NCornerNodes <<
330  " NSides " << NSides;
331  LOGPZ_DEBUG(logger,sout.str());
332  }
333 
334  if (logger->isDebugEnabled())
335  {
336  std::stringstream sout;
337  sout << "Testing smallsides\n";
338  TPZStack<int> smallsides;
339  int is;
340  for(is=0; is<NSides; is++)
341  {
342  LowerDimensionSides(is,smallsides);
343  sout << " side " << is << " allsmallsides " << smallsides << endl;
344  int id;
345  for(id=0; id < SideDimension(is); id++)
346  {
347  smallsides.Resize(0);
348  LowerDimensionSides(is,smallsides,id);
349  sout << " side " << is << " targetdimension " << id << " smallsides " << smallsides << endl;
350  }
351  }
352  LOGPZ_DEBUG(logger,sout.str());
353  }
354 
355  // static void HigherDimensionSides(int side, TPZStack<int> &high);
356  if (logger->isDebugEnabled())
357  {
358  std::stringstream sout;
359  sout << "Testing HigherDimensionSides\n";
360  int is ;
361  for(is=0; is<NSides; is++)
362  {
363  TPZStack<int> high;
364  HigherDimensionSides(is,high);
365  sout << "Side " << is << " higherdimensionsides " << high << endl;
366  }
367  LOGPZ_DEBUG(logger,sout.str());
368  }
369  // static int NSideNodes(int side);
370  if (logger->isDebugEnabled())
371  {
372  std::stringstream sout;
373  sout << "Testing NSideNodes\n";
374  int is;
375  for(is=0; is<NSides; is++)
376  {
377  sout << "Side " << is << " NSideNodes " << NSideNodes(is) << endl;
378  }
379  LOGPZ_DEBUG(logger,sout.str());
380  }
384  // static int SideNodeLocId(int side, int node);
385  if (logger->isDebugEnabled())
386  {
387  std::stringstream sout;
388  sout << "Testing SideNodeLocId\n";
389  int is;
390  for(is=0; is<NSides; is++)
391  {
392  int nsnodes = NSideNodes(is);
393  int sn;
394  sout << "Side " << is << " sidenodelocid ";
395  for(sn=0; sn<nsnodes; sn++)
396  {
397  sout << SideNodeLocId(is,sn) << " ";
398  }
399  sout << endl;
400  }
401  LOGPZ_DEBUG(logger,sout.str());
402  }
403  // static void CenterPoint(int side, TPZVec<REAL> &center);
404  if (logger->isDebugEnabled())
405  {
406  std::stringstream sout;
407  sout << "Testing CenterPoint\n";
408  int is;
409  for(is=0; is<NSides; is++)
410  {
411  TPZManVector<REAL> center(Dimension);
412  CenterPoint(is,center);
413  sout << "Side " << is << center << endl;
414  }
415  LOGPZ_DEBUG(logger,sout.str());
416  }
417  // static REAL RefElVolume(){return 2.0*TFather::RefElVolume();}
418  if (logger->isDebugEnabled())
419  {
420  std::stringstream sout;
421  sout << "Testing RefelVolume\n";
422  sout << "Volume " << RefElVolume() << endl;
423  LOGPZ_DEBUG(logger,sout.str());
424  }
425  // static int SideDimension(int side);
426  if (logger->isDebugEnabled())
427  {
428  std::stringstream sout;
429  sout << "Testing SideDimension\n";
430  int is;
431  for(is=0; is<NSides; is++)
432  {
433  sout << "Side " << is << " sidedimension " << SideDimension(is) << endl;
434  }
435  LOGPZ_DEBUG(logger,sout.str());
436  }
437 
438  // static TPZTransform<> SideToSideTransform(int sidefrom, int sideto);
439  if (logger->isDebugEnabled())
440  {
441  std::stringstream sout;
442  sout << "Testing SideToSideTransform\n";
443  int is;
444  for(is=0; is<NSides; is++)
445  {
446  TPZStack<int> lower;
447  LowerDimensionSides(is,lower);
448  int nl = lower.NElements();
449  int il;
450  for(il=0; il<nl; il++)
451  {
452  TPZTransform<> tr;
453  tr = SideToSideTransform(is,lower[il]);
454  sout << "Transform from " << is << " to side " << lower[il] << " " << tr;
455  }
456  }
457  LOGPZ_DEBUG(logger,sout.str());
458  }
459  // static TPZTransform<> TransformElementToSide(int side);
460  if (logger->isDebugEnabled())
461  {
462  std::stringstream sout;
463  sout << "Testing TransformElementToSide\n";
464  int is;
465  for(is=0; is<NSides; is++)
466  {
467  TPZTransform<> tr = TransformElementToSide(is);
468  sout << "side " << is << " transform " << tr;
469  }
470  LOGPZ_DEBUG(logger,sout.str());
471  }
472  // static TPZTransform<> TransformSideToElement(int side);
473  if (logger->isDebugEnabled())
474  {
475  std::stringstream sout;
476  sout << "Testing TransformSideToElement\n";
477  int is;
478  for(is=0; is<NSides; is++)
479  {
480  sout << "side " << is << " transform " << TransformSideToElement(is);
481  }
482  LOGPZ_DEBUG(logger,sout.str());
483  }
484 
485  // static TPZIntPoints *CreateSideIntegrationRule(int side, int order);
486  if (logger->isDebugEnabled())
487  {
488  std::stringstream sout;
489  sout << "Testing CreateSideIntegrationRule order 3\n";
490  int is;
491  for(is=0; is<NSides; is++)
492  {
493  TPZIntPoints *integ = CreateSideIntegrationRule(is,3);
494  sout << "side " << is << endl;
495  integ->Print(sout);
496  delete integ;
497  }
498  LOGPZ_DEBUG(logger,sout.str());
499  }
500 
501  // static std::string StrType() ;//{ return EOned;}
502  if (logger->isDebugEnabled())
503  {
504  std::stringstream sout;
505  sout << "Testing StrType\n";
506  sout << StrType();
507  LOGPZ_DEBUG(logger,sout.str());
508  }
509 
510  // static std::string StrType(int side);
511  if (logger->isDebugEnabled())
512  {
513  std::stringstream sout;
514  sout << "Testing StrType(side)\n";
515  int is;
516  for(is=0; is<NSides; is++)
517  {
518  sout << "side " << is << " type " << StrType(is) << endl;
519  }
520  LOGPZ_DEBUG(logger,sout.str());
521  }
522 
523  // static void MapToSide(int side, TPZVec<REAL> &InternalPar, TPZVec<REAL> &SidePar, TPZFMatrix<REAL> &JacToSide);
524  if (logger->isDebugEnabled())
525  {
526  std::stringstream sout;
527  TPZManVector<REAL> par(Dimension,sqrt(2.));
528  sout << "Testing MapToSide internalpar = " << par << "\n";
529  int is;
530  for(is=0; is<NSides; is++)
531  {
532  TPZManVector<REAL> sidepar(SideDimension(is));
533  TPZFMatrix<REAL> jactoside (SideDimension(is),Dimension);
534  MapToSide(is,par,sidepar,jactoside);
535  sout << "side " << is << " sidepar " << sidepar << " jactoside " << jactoside << endl;
536  }
537  LOGPZ_DEBUG(logger,sout.str());
538  }
539 
540  // static int NSides;
541  if (logger->isDebugEnabled())
542  {
543  std::stringstream sout;
544  sout << "Testing NumSides\n";
545  sout << "nsides " << NumSides();
546  LOGPZ_DEBUG(logger,sout.str());
547  }
548 
549  // static int NContainedSides(int side);
550  if (logger->isDebugEnabled())
551  {
552  std::stringstream sout;
553  sout << "Testing NContainedSides(side)\n";
554  int is;
555  for(is=0; is<NSides; is++)
556  {
557  sout << "side " << is << " nsideconnects " << NContainedSides(is) << endl;
558  }
559  LOGPZ_DEBUG(logger,sout.str());
560  }
561  // static int ContainedSideLocId(int side, int c);
562  if (logger->isDebugEnabled())
563  {
564  std::stringstream sout;
565  sout << "Testing ContainedSideLocId(is,ic)\n";
566  int is,ic;
567  for(is=0; is<NSides; is++)
568  {
569  int nc = NContainedSides(is);
570  sout << "side " << is << " connects ";
571  for(ic=0; ic<nc; ic++)
572  {
573  sout << ContainedSideLocId(is,ic) << " ";
574  }
575  sout << endl;
576  }
577 
578  LOGPZ_DEBUG(logger,sout.str());
579  }
580 
581 
582 #endif
583  }
584 
585 
586  template<>
587  REAL Pr< pztopology::TPZPoint >::RefElVolume()
588  {
589  return 2.0L;
590  }
591 
592  template class Pr<TPZPoint>;
593  template class Pr<Pr<TPZPoint> >;
594 
595 }
596 
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZTransform::Mult
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
tpzpoint.h
Contains the TPZPoint class which defines the topology of a point.
MElementType_Name
std::string MElementType_Name(MElementType elType)
Returns the name of the element type.
Definition: pzeltype.h:127
TPZIntPoints::PrismExtend
virtual TPZIntPoints * PrismExtend(int order)=0
pzrefine::fatherside
static int fatherside[8][21]
Definition: pzrefprism.cpp:303
PrismExtend.h
Contains the Pr class(template) which defines the Prismatic extension of a topology.
TPZIntPoints::Print
virtual void Print(std::ostream &out) const
Prints information of the cubature rule.
Definition: pzquad.cpp:38
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
TPZIntPoints
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
pztopology::highsides
static int highsides[27][7]
For each side was stored the sides connected with it but of the higher dimension. ...
Definition: tpzcube.cpp:110
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
pztopology
Groups all classes defining the structure of the master element.
Definition: PrismExtend.cpp:15
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZTransform::Sum
const TPZFMatrix< T > & Sum() const
Definition: pztrnsform.h:66
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZTransform::Multiply
TPZTransform< T > Multiply(TPZTransform< T > &right)
Multiply the transformation object (to the right) with right (Multiplying matrices) ...
Definition: pztrnsform.cpp:112
pztopology::Pr
Defines the Prismatic extension of a topology. Topology.
Definition: PrismExtend.h:25
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZTransform
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
TPZTransform::Apply
void Apply(TPZVec< T > &vectorin, TPZVec< T > &vectorout)
Transforms the vector.
Definition: pztrnsform.cpp:121
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
Generated on Sun Aug 16 2020 11:17:11 for NeoPZ by   doxygen 1.8.13