NeoPZ
tpzline.cpp
Go to the documentation of this file.
1 
6 #include "tpzline.h"
7 
8 #include "pzerror.h"
9 #include "pzreal.h"
10 #include "pztrnsform.h"
11 #include "pzquad.h"
12 #include "pzeltype.h"
13 
14 #ifdef _AUTODIFF
15 #include "fad.h"
16 #endif
17 
18 
19 #include "pzlog.h"
20 
21 #ifdef LOG4CXX
22 static LoggerPtr logger(Logger::getLogger("pz.topology.pzline"));
23 #endif
24 
25 using namespace std;
26 
27 namespace pztopology {
28 
29  static int nhighdimsides[3] = {1,1,0};
30 
31  static int sidedimension[3] = {0,0,1};
32 
33  static int highsides[3][1] = {
34  {2},
35  {2},
36  {0}
37  };
38 
39  static REAL sidetosidetransforms[3][1][4][3] = {
40  {
41  {{0,0,0},{0,0,0},{0,0,0},{-1,0,0}}
42  },
43  {
44  {{0,0,0},{0,0,0},{0,0,0},{1,0,0}}
45  },
46  {
47  {{0,0,0},{0,0,0},{0,0,0},{0,0,0}}}
48  };
49 
50  static REAL MidSideNode[3][1] = {{-1.},{1.},{0.}};
51 
52  static int nsidenodes[3] = {1,1,2};
53 
54  int TPZLine::fPermutations [2][3] =
55  {
56  {0,1,2},
57  {1,0,2}
58  };
59 
60  int TPZLine::NBilinearSides()
61  {return 0;}
62 
63  static int vectorsideorder [3] = {0,1,2};
64 
65  static int bilinearounao [3] = {0,0,1};
66  static int direcaoksioueta [3] = {0,0,0};
67 
68  template<class T>
69  inline void TPZLine::TShape(const TPZVec<T> &loc,TPZFMatrix<T> &phi,TPZFMatrix<T> &dphi) {
70  T x = loc[0];
71  phi(0,0) = (1.0-x)/2.;
72  phi(1,0) = (1.0+x)/2.;
73  dphi(0,0) = -0.5;
74  dphi(0,1) = 0.5;
75  }
76 
77  template<class T>
78  void TPZLine::BlendFactorForSide(const int &side, const TPZVec<T> &xi, T &blendFactor,
79  TPZVec<T> &blendFactorDxi){
80  const REAL tol = pztopology::GetTolerance();
81 #ifdef PZDEBUG
82  std::ostringstream sout;
83  if(side < NCornerNodes || side >= NSides){
84  sout<<"The side\t"<<side<<"is invalid. Aborting..."<<std::endl;
85  PZError<<std::endl<<sout.str()<<std::endl;
86  DebugStop();
87  }
88 
89  if(!TPZLine::IsInParametricDomain(xi,tol)){
90  sout<<"The method BlendFactorForSide expects the point xi to correspond to coordinates of a point";
91  sout<<" inside the parametric domain. Aborting...";
92  PZError<<std::endl<<sout.str()<<std::endl;
93  #ifdef LOG4CXX
94  LOGPZ_FATAL(logger,sout.str().c_str());
95  #endif
96  DebugStop();
97  }
98 #endif
99  TPZFNMatrix<4,T> phi(NCornerNodes,1);
100  TPZFNMatrix<8,T> dphi(Dimension,NCornerNodes);
101  TPZLine::TShape(xi,phi,dphi);
102  blendFactorDxi.Resize(TPZLine::Dimension, (T) 0);
103  switch(side){
104  case 0:
105  blendFactor = phi(0,0);
106  blendFactorDxi[0] = dphi(0,0);
107  break;
108  case 1:
109  blendFactor = phi(1,0);
110  blendFactorDxi[0] = dphi(0,1);
111  break;
112  case 2:
113  blendFactor = 1;
114  blendFactorDxi[0] = 0;
115  break;
116  }
117  return;
118  }
119 
120  int TPZLine::NSideNodes(int side)
121  {
122  return nsidenodes[side];
123  }
124 
125  int TPZLine::NumSides(int dimension) {
126  if(dimension<0 || dimension> 1) {
127  PZError << "TPZLine::NumSides. Bad parameter i.\n";
128  return 0;
129  }
130  if(dimension==0) return 2;
131  if(dimension==1) return 1;
132  return -1;
133  }
134 
135 
136  void TPZLine::LowerDimensionSides(int side,TPZStack<int> &smallsides)
137  {
138  smallsides.Resize(0);
139  int nsidecon = NContainedSides(side);
140  int is;
141  for(is=0; is<nsidecon-1; is++)
142  smallsides.Push(ContainedSideLocId(side,is));
143  }
144 
145  void TPZLine::LowerDimensionSides(int side,TPZStack<int> &smallsides, int DimTarget)
146  {
147  smallsides.Resize(0);
148  int nsidecon = NContainedSides(side);
149  for(int is = 0; is < nsidecon - 1; is++) {
150  if (SideDimension(ContainedSideLocId(side,is)) == DimTarget) smallsides.Push(ContainedSideLocId(side,is));
151  }
152  }
153 
154  bool TPZLine::IsInParametricDomain(const TPZVec<REAL> &pt, REAL tol) {
155  const REAL qsi = pt[0];
156  if( fabs(qsi) <= 1. + tol){
157  return true;
158  }
159  else{
160  return false;
161  }
162  }//method
163 
164  template<class T>
165  bool TPZLine::CheckProjectionForSingularity(const int &side, const TPZVec<T> &xiInterior) {
166  return true;
167  }
168 
169  template<class T>
170  void TPZLine::MapToSide(int side, TPZVec<T> &InternalPar, TPZVec<T> &SidePar, TPZFMatrix<T> &JacToSide) {
171  TPZTransform<> TransfR = SideToSideTransform(NSides - 1, side);
172  TPZTransform<T> Transf;
173  Transf.CopyFrom(TransfR);
174  SidePar.Resize(SideDimension(side));
175  Transf.Apply(InternalPar,SidePar);
176 
177  JacToSide = Transf.Mult();
178  }
179 
180  void TPZLine::ParametricDomainNodeCoord(int node, TPZVec<REAL> &nodeCoord)
181  {
182  if(node > NCornerNodes)
183  {
184  DebugStop();
185  }
186  nodeCoord.Resize(Dimension, 0.);
187  switch (node) {
188  case (0):
189  {
190  nodeCoord[0] = -1.;
191  break;
192  }
193  case (1):
194  {
195  nodeCoord[0] = 1.;
196  break;
197  }
198  default:
199  {
200  DebugStop();
201  break;
202  }
203  }
204  }
205 
206  void TPZLine::HigherDimensionSides(int side, TPZStack<int> &high)
207  {
208  if(side <0 || side >= NSides) {
209  PZError << "TPZLine::HigherDimensionSides side "<< side << endl;
210  }
211  int is;
212  for(is=0; is<nhighdimsides[side]; is++) high.Push(highsides[side][is]);
213 
214  }
215 
216  int TPZLine::SideNodeLocId(int side, int node)
217  {
218  if(side <2 && node == 0) return side;
219  if(side == 2 && node <2) return node;
220  PZError << "TPZLine::SideNodeLocId inconsistent side or node " << side
221  << ' ' << node << endl;
222  return -1;
223  }
224  void TPZLine::CenterPoint(int side, TPZVec<REAL> &center) {
225  if (center.size()!=Dimension) {
226  DebugStop();
227  }
228  int i;
229  for(i=0; i<Dimension; i++) {
230  center[i] = MidSideNode[side][i];
231  }
232  }
233 
234 
236  void TPZLine::RandomPoint(TPZVec<REAL> &pt)
237  {
238  for(int i=0; i<1; i++)
239  {
240  REAL val = -1. + 2.*(REAL) rand() / (RAND_MAX);
241  pt[i] = val;
242  }
243  }
244 
245  int TPZLine::SideDimension(int side) {
246  if(side<0 || side >= NSides) {
247  PZError << "TPZLine::SideDimension side " << side << endl;
248  return -1;
249  }
250  return sidedimension[side];
251  }
252 
253  TPZTransform<> TPZLine::SideToSideTransform(int sidefrom, int sideto)
254  {
255  if(sidefrom <0 || sidefrom >= NSides || sideto <0 || sideto >= NSides) {
256  PZError << "TPZLine::HigherDimensionSides sidefrom "<< sidefrom <<
257  ' ' << sideto << endl;
258  return TPZTransform<>(0);
259  }
260  if(sidefrom == sideto) {
261  return TPZTransform<>(sidedimension[sidefrom]);
262  }
263  if(sidefrom == NSides-1) {
264  return TransformElementToSide(sideto);
265  }
266 
267  if (sideto == NSides -1) {
268  return TransformSideToElement(sidefrom);
269  }
270 
271  int nhigh = nhighdimsides[sidefrom];
272  int is;
273  for(is=0; is<nhigh; is++) {
274  if(highsides[sidefrom][is] == sideto) {
275  int dfr = sidedimension[sidefrom];
276  int dto = sidedimension[sideto];
277  TPZTransform<> trans(dto,dfr);
278  int i,j;
279  for(i=0; i<dto; i++) {
280  for(j=0; j<dfr; j++) {
281  trans.Mult()(i,j) = sidetosidetransforms[sidefrom][is][j][i];
282  }
283  trans.Sum()(i,0) = sidetosidetransforms[sidefrom][is][3][i];
284  }
285  return trans;
286  }
287  }
288  PZError << "TPZLine::SideToSideTransform highside not found sidefrom "
289  << sidefrom << ' ' << sideto << endl;
290  return TPZTransform<>(0);
291  }
292 
293  TPZTransform<> TPZLine::TransformElementToSide(int side){
294 
295  if(side<0 || side>2){
296  PZError << "TPZLine::TransformElementToSide called with side error\n";
297  return TPZTransform<>(0,0);
298  }
299 
300  // int sidedim = SideDimension(side);
301  TPZTransform<> t(SideDimension(side),1);//t(dimto,2)
302  t.Mult().Zero(); //TPZGeoElQ2d *gq;
303  t.Sum().Zero();//int dimto = gq->SideDimension(side);
304 
305  switch(side){
306  case 0:
307  case 1:
308  return t;
309  case 2:
310  t.Mult()(0,0) = 1.0;//par. var.
311  return t;
312  }
313  return TPZTransform<>(0,0);
314  }
315 
316  TPZTransform<> TPZLine::TransformSideToElement(int side){
317 
318  if(side<0 || side>2){
319  PZError << "TPZLine::TransformSideToElement side out range\n";
320  return TPZTransform<>(0,0);
321  }
322  int sidedim = 1;
323  if(side <2) sidedim = 0;
324 
325  TPZTransform<> t(1,sidedim);
326  t.Mult().Zero();
327  t.Sum().Zero();
328 
329  switch(side){
330  case 0:
331  t.Sum()(0,0) = -1.0;
332  return t;
333  case 1:
334  t.Sum()(0,0) = 1.0;
335  return t;
336  case 2:
337  t.Mult()(0,0) = 1.0;
338  return t;
339  }
340  return TPZTransform<>(0,0);
341  }
342 
343 
344  TPZIntPoints *TPZLine::CreateSideIntegrationRule(int side, int order) {
345 
346  if(side<0 || side>2) {
347  PZError << "TPZLine::CreateSideIntegrationRule wrong side " << side << endl;
348  return 0;
349  }
350  if(side != 2) return new TPZInt1Point(order); // sides 0 and 1 are vertices (corners)
351  return new IntruleType(order);
352 
353 
354  }
355 
356 
357  MElementType TPZLine::Type()
358  {
359  return EOned;
360  }
361 
362  MElementType TPZLine::Type(int side)
363  {
364  switch(side) {
365  case 0:
366  case 1:
367  return EPoint;
368  case 2:
369  return EOned;
370  default:
371  return ENoType;
372  }
373  }
374 
375 
376  int TPZLine::NumSides() {
377  return 3;
378  }
379 
380 
381  int TPZLine::NContainedSides(int side) {
382  if(side==0 || side==1) return 1;
383  else if(side==2) return 3;
384  PZError << "TPZLine::NContainedSides. Bad parameter side = " << side << " .\n";
385  DebugStop();
386  return 0;
387  }
388 
389  int TPZLine::ContainedSideLocId(int side,int c) {
390  switch(side) {
391  case 0:
392  case 1:
393  if(c != 0)
394  {
395  PZError << "TPZLine::ContainedSideLocId, connect = " << c << endl;
396  }
397  return side;
398  case 2:
399  return c;
400  default:
401  PZError << "TPZLine::ContainedSideLocId called with side = " << side << endl;
402  return 0;
403  }
404  }
405 
412  int TPZLine::GetTransformId(TPZVec<int64_t> &id)
413  {
414  return id[0] < id[1] ? 0 : 1;
415  }
416 
423  int TPZLine::GetTransformId(int side, TPZVec<int64_t> &id)
424  {
425  switch (side) {
426  case 0:
427  case 1:
428  return 0;
429  break;
430  case 2:
431  return id[0] < id[1] ? 0 : 1;
432  default:
433  break;
434  }
435  LOGPZ_ERROR(logger,"Wrong input parameter")
436  return -1;
437  }
438 
439 
447  void TPZLine::GetSideHDivPermutation(int transformationid, TPZVec<int> &permgather)
448  {
449  permgather.Resize(3);
450  for (int i=0; i<3; i++)
451  {
452  permgather[i] = fPermutations[transformationid][i];
453  }
454  /*
455  switch (side) {
456  case 0:
457  case 1:
458  permgather[0] = 0;
459  break;
460  case 2:
461  if(id[0]<id[1])
462  {
463  permgather[0] = 0;
464  permgather[1] = 1;
465  permgather[2] = 2;
466  }
467  else
468  {
469  permgather[0] = 1;
470  permgather[1] = 0;
471  permgather[2] = 2;
472  }
473 
474  default:
475  break;
476  }
477  LOGPZ_ERROR(logger,"Wrong input parameter")*/
478  }
479 
480 
481  void TPZLine::ComputeDirections(int side, TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions, TPZVec<int> &sidevectors)
482  {
483 
484  if(gradx.Cols()!=1)
485  {
486  DebugStop();
487  }
488  TPZFMatrix<REAL> bvec = gradx;
489  REAL norma = 0.0;
490  for (int i = 0; i<3; i++)
491  {
492  norma += bvec(i,0)*bvec(i,0);
493  }
494  norma = sqrt(norma);
495 
496  directions.Redim(3, 1);
497  for (int i=0; i<3; i++)
498  {
499  directions(i,0) = bvec(i,0)/norma;
500  }
501  sidevectors.Resize(1);
502  sidevectors[0] = -1; // ??????????????????????
503 
504  }
505 
506  void TPZLine::GetSideHDivDirections(TPZVec<int> &sides, TPZVec<int> &dir, TPZVec<int> &bilounao)
507  {
508  int nsides = NumSides();
509 
510  sides.Resize(nsides);
511  dir.Resize(nsides);
512  bilounao.Resize(nsides);
513 
514 
515  for (int is = 0; is<nsides; is++)
516  {
517  sides[is] = vectorsideorder[is];
518  dir[is] = direcaoksioueta[is];
519  bilounao[is] = bilinearounao[is];
520  }
521  }
522 
523  void TPZLine::GetSideHDivDirections(TPZVec<int> &sides, TPZVec<int> &dir, TPZVec<int> &bilounao, TPZVec<int> &sidevectors)
524  {
525  int nsides = NumSides();
526 
527  sides.Resize(nsides);
528  dir.Resize(nsides);
529  bilounao.Resize(nsides);
530 
531  for (int is = 0; is<nsides; is++)
532  {
533  sides[is] = vectorsideorder[is];
534  dir[is] = direcaoksioueta[is];
535  bilounao[is] = bilinearounao[is];
536  sidevectors[is] = is;
537  }
538  }
539 
540  int TPZLine::ClassId() const{
541  return Hash("TPZLine");
542  }
543 
544  void TPZLine::Read(TPZStream& buf, void* context) {
545 
546  }
547 
548  void TPZLine::Write(TPZStream& buf, int withclassid) const {
549 
550  }
551 
552 }
553 
554 /**********************************************************************************************************************
555  * The following are explicit instantiation of member function template of this class, both with class T=REAL and its
556  * respective FAD<REAL> version. In other to avoid potential errors, always declare the instantiation in the same order
557  * in BOTH cases. @orlandini
558  **********************************************************************************************************************/
559 template bool pztopology::TPZLine::CheckProjectionForSingularity<REAL>(const int &side, const TPZVec<REAL> &xiInterior);
560 
561 template void pztopology::TPZLine::MapToSide<REAL>(int side, TPZVec<REAL> &InternalPar, TPZVec<REAL> &SidePar, TPZFMatrix<REAL> &JacToSide);
562 
563 template void pztopology::TPZLine::BlendFactorForSide<REAL>(const int &, const TPZVec<REAL> &, REAL &, TPZVec<REAL> &);
564 
565 template void pztopology::TPZLine::TShape<REAL>(const TPZVec<REAL> &loc,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi);
566 
567 template void pztopology::TPZLine::ComputeHDivDirections<REAL>(TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions);
568 #ifdef _AUTODIFF
569 
570 template bool pztopology::TPZLine::CheckProjectionForSingularity<Fad<REAL>>(const int &side, const TPZVec<Fad<REAL>> &xiInterior);
571 
572 template void pztopology::TPZLine::MapToSide<Fad<REAL> >(int side, TPZVec<Fad<REAL> > &InternalPar, TPZVec<Fad<REAL> > &SidePar, TPZFMatrix<Fad<REAL> > &JacToSide);
573 
574 template void pztopology::TPZLine::BlendFactorForSide<Fad<REAL>>(const int &, const TPZVec<Fad<REAL>> &, Fad<REAL> &,
575  TPZVec<Fad<REAL>> &);
576 template void pztopology::TPZLine::TShape<Fad<REAL>>(const TPZVec<Fad<REAL>> &loc,TPZFMatrix<Fad<REAL>> &phi,TPZFMatrix<Fad<REAL>> &dphi);
577 
578 template void pztopology::TPZLine::ComputeHDivDirections<Fad<REAL>>(TPZFMatrix<Fad<REAL>> &gradx, TPZFMatrix<Fad<REAL>> &directions);
579 #endif
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
pzquad.h
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
pztopology::bilinearounao
static int bilinearounao[3]
Definition: tpzline.cpp:65
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.
dimension
clarg::argInt dimension("-d", "Matrices dimension M x M", 1000)
ENoType
Definition: pzeltype.h:71
TPZTransform::CopyFrom
void CopyFrom(const TPZTransform< REAL > &cp)
Create a copy form a real transformation.
Definition: pztrnsform.cpp:62
TPZTransform::Mult
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
pzerror.h
Defines PZError.
Fad
Definition: fad.h:54
pztopology::MidSideNode
static REAL MidSideNode[3][1]
Definition: tpzline.cpp:50
pzeltype.h
Defines enum MElementType and contains the implementation of MElementType_NNodes(...) functions.
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
TPZInt1d
Handles the numerical integration for one-dimensional problems. Numerical Integration.
Definition: pzquad.h:36
LOGPZ_FATAL
#define LOGPZ_FATAL(A, B)
Define log for fatal errors (cout)
Definition: pzlog.h:95
TPZIntPoints
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
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
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
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
pztopology::vectorsideorder
static int vectorsideorder[3]
Definition: tpzline.cpp:63
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
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
pztopology::nhighdimsides
static int nhighdimsides[3]
Definition: tpzline.cpp:29
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
pztopology::nsidenodes
static int nsidenodes[3]
Definition: tpzline.cpp:52
EPoint
Definition: pzeltype.h:54
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
pztopology::highsides
static int highsides[3][1]
Definition: tpzline.cpp:33
TPZInt1Point
Integration rule for one point. Numerical Integration.
Definition: pzquad.h:478
MElementType
MElementType
Define the element types.
Definition: pzeltype.h:52
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
pztrnsform.h
Contains the TPZTransform<> class which implements an affine transformation between points in paramet...
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
pztopology::sidetosidetransforms
static REAL sidetosidetransforms[3][1][4][3]
Definition: tpzline.cpp:39
pzreal.h
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
TPZTransform
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
tpzline.h
Contains the TPZLine class which defines the topology of a line element.
EOned
Definition: pzeltype.h:55
pztopology::sidedimension
static int sidedimension[3]
Definition: tpzline.cpp:31
TPZTransform::Apply
void Apply(TPZVec< T > &vectorin, TPZVec< T > &vectorout)
Transforms the vector.
Definition: pztrnsform.cpp:121
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
pztopology::direcaoksioueta
static int direcaoksioueta[3]
Definition: tpzline.cpp:66
Generated on Sun Aug 16 2020 11:17:11 for NeoPZ by   doxygen 1.8.13