NeoPZ
pzshapelinear.cpp
Go to the documentation of this file.
1 
6 #include "pzshapelinear.h"
7 #include "pzshapepoint.h"
8 #include "pzerror.h"
9 #include "pzreal.h"
10 using namespace std;
11 
12 namespace pzshape {
13 
14 
15  void TPZShapeLinear::Chebyshev(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi){
16  // Quadratic or higher shape functions
17  if(num <= 0) return;
18  phi.Zero();
19  dphi.Zero();
20  phi.Put(0,0,1.0);
21  dphi.Put(0,0, 0.0);
22  if(num == 1) return;
23  phi.Put(1,0, x);
24  dphi.Put(0,1, 1.0);
25  int ord;
26  // dphi.Print("DphisAntes = ",std::cout,EMathematicaInput);
27 
28  for(ord = 2;ord<num;ord++) {
29  phi.Put(ord,0, 2.0*x*phi(ord-1,0) - phi(ord-2,0));
30  dphi.Put(0,ord, 2.0*x*dphi(0,ord-1) + 2.0*phi(ord-1,0) - dphi(0,ord-2));
31  }
32  // dphi.Print("DphisDepois = ",std::cout,EMathematicaInput);
33  }
34 
35  void TPZShapeLinear::Expo(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi){
36  // Quadratic or higher shape functions
37  if(num <= 0) return;
38  phi.Put(0,0,1.0);
39  dphi.Put(0,0, 0.0);
40  if(num == 1) return;
41  phi.Put(1,0, x);
42  dphi.Put(0,1, 1.0);
43  int ord;
44  for(ord = 2;ord<num;ord++) {
45  phi.Put(ord,0, x*phi(ord-1,0));
46  dphi.Put(0,ord, ord*phi(ord-1,0));
47  }
48  }
49 
50  void TPZShapeLinear::Legendre(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi){
51 
52  // Quadratic or higher shape functions
53  if (num <= 0) return;
54  phi.Put(0, 0, 1.0);
55  dphi.Put(0, 0, 0.0);
56 
57  if (num == 1) return;
58 
59  phi.Put(1, 0, x);
60  dphi.Put(0, 1, 1.0);
61 
62  int ord;
63  //Aqui fica diferente do Chebyshev
64  REAL ord_real, value;
65  for (ord = 2; ord < num; ord++)
66  {
67  //casting int ord to REAL ord_real
68  ord_real = (REAL)ord;
69  //computing the ord_th function
70  value = ( (2.0 * (ord_real - 1.0) + 1.0) * x * phi(ord - 1, 0) - (ord_real - 1.0) * phi(ord - 2 , 0) ) / (ord_real);
71  phi.Put(ord, 0, value);
72 
73  //computing the ord_th function's derivative
74  value = (2.0 * (ord_real - 1.0) + 1.0) * phi(ord - 1, 0) + dphi(0, ord - 2);
75  dphi.Put(0, ord, value);
76  }
77 
78 #ifdef PZDEBUG
79  int printing = 0;
80  if (printing){
81  cout << "Legendre" << endl;
82  for(ord = 0; ord < num; ord++)
83  {
84  cout << "x = " << x << endl;
85  cout << "phi(" << ord << ", 0) = " << phi(ord, 0) << endl;
86  cout << "dphi(0, " << ord << " = " << dphi(0, ord) << endl;
87  cout << endl;
88  }
89  }
90 #endif
91 
92  } //end of method
93 
94  void TPZShapeLinear::Legendre(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi, int nderiv){
95 
96  // Quadratic or higher shape functions
97  if (num <= 0) return;
98  phi.Put(0, 0, 1.0);
99  dphi.Put(0, 0, 0.0);
100 
101  int ideriv;
102  for (ideriv = 1; ideriv < nderiv; ideriv++)
103  dphi.Put(ideriv, 0, 0.0);
104 
105 
106  if (num == 1) return;
107 
108  phi.Put(1, 0, x);
109  dphi.Put(0, 1, 1.0);
110 
111  for (ideriv = 1; ideriv < nderiv; ideriv++)
112  dphi.Put(ideriv, 1, 0.0);
113 
114  int ord;
115  //Aqui fica diferente do Chebyshev
116  REAL ord_real, value;
117  for (ord = 2; ord < num; ord++)
118  {
119  //casting int ord to REAL ord_real
120  ord_real = (REAL)ord;
121  //computing the ord_th function
122  value = ( (2.0 * (ord_real - 1.0) + 1.0) * x * phi(ord - 1, 0) - (ord_real - 1.0) * phi(ord - 2 , 0) ) / (ord_real);
123  phi.Put(ord, 0, value);
124 
125  //computing the ord_th function's derivative
126  value = (2.0 * (ord_real - 1.0) + 1.0) * phi(ord - 1, 0) + dphi(0, ord - 2);
127  dphi.Put(0, ord, value);
128 
129  for (ideriv = 1; ideriv < nderiv; ideriv++){
130  value = (2.0 * (ord_real - 1.0) + 1.0) * dphi(ideriv - 1, ord - 1) + dphi(ideriv, ord - 2);
131  dphi.Put(ideriv, ord, value);
132  }
133 
134  }
135 
136 #ifdef PZDEBUG
137  int printing = 0;
138  if (printing){
139  cout << "Legendre" << endl;
140  for(ord = 0; ord < num; ord++)
141  {
142  cout << "x = " << x << endl;
143  cout << "phi(" << ord << ", 0) = " << phi(ord, 0) << endl;
144  cout << "dphi(0, " << ord << " = " << dphi(0, ord) << endl;
145  cout << endl;
146  }
147  }
148 #endif
149 
150  } //end of method
151 
155  REAL TPZShapeLinear::JacobiA = 0.5;
156  REAL TPZShapeLinear::JacobiB = 0.5;
157 
158  void TPZShapeLinear::Jacobi(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
159  // Quadratic or higher shape functions
160  REAL a = JacobiA, b = JacobiB;
161  if(num <= 0) return;
162  phi.Put(0,0,1.0);
163  dphi.Put(0,0, 0.0);
164  if(num == 1) return;
165  phi.Put(1,0, 0.5 * (a - b + (2.0 + a + b) * x));
166  dphi.Put(0,1, 0.5*(num + a + b + 1)*phi(num-1,0) );
167  // Following http://mathworld.wolfram.com/JacobiPolynomial.html
168  REAL A_ab=0.0, B_ab=0.0, C_ab=0.0;
169  for(int ord = 2;ord<num;ord++) {
170  //Computing Coefficients for three terms recursion relation
171  A_ab = ((2.0*ord+a+b+1)*(2.0*ord+a+b+2))/(2.0*(ord+1)*(ord+a+b+1));
172  B_ab = ((b*b-a*a)*(2.0*ord+a+b+1))/(2.0*(ord+1)*(ord+a+b+1)*(2.0*ord+a+b));
173  C_ab = ((ord+a)*(ord+b)*(2.0*ord+a+b+2))/((ord+1)*(ord+a+b+1)*(2.0*ord+a+b));
174  phi.Put(ord,0, (A_ab*x - B_ab)*phi(ord-1,0) - C_ab*phi(ord-2,0));
175  dphi.Put(0,ord, 0.5*(ord + a + b + 1)*phi(ord-1,0));
176  }
177  } //end of method
178 
179  void TPZShapeLinear::Hermite(REAL x,int num,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
180  // Quadratic or higher shape functions (for physicists)
181  if(num <= 0) return;
182  phi.Put(0,0,1.0);
183  dphi.Put(0,0, 0.0);
184  if(num == 1) return;
185  phi.Put(1,0, 2.0*x);
186  dphi.Put(0,1, 2.0);
187  // Following http://mathworld.wolfram.com/HermitePolynomial.html
188  for(int ord = 2;ord<num;ord++) {
189  phi.Put(ord,0, 2.0*x*phi(ord-1,0) - (2.0*(REAL)ord-1.0)*phi(ord-2,0));
190  dphi.Put(0,ord, (2.0*(REAL)ord-1.0)*phi(ord-2,0));
191  }
192  } //end of method
193 
194  // Setting Chebyshev polynomials as orthogonal sequence generating shape functions
195  void (*TPZShapeLinear::fOrthogonal)(REAL, int, TPZFMatrix<REAL> &, TPZFMatrix<REAL> &) = TPZShapeLinear::Chebyshev;
196 
203  void TPZShapeLinear::ShapeGenerating(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi)
204  {
205 
206  phi(2,0) = phi(0,0)*phi(1,0);
207  dphi(0,2) = dphi(0,0)*phi(1,0)+phi(0,0)*dphi(0,1);
208 
209  phi(2,0) *= 4.;
210  dphi(0,2) *= 4.;
211 
212  }
213 
214  void TPZShapeLinear::Shape(TPZVec<REAL> &x,TPZVec<int64_t> &id, TPZVec<int> &order,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
215  // num = number of functions to compute
216 #ifndef NODEBUG
217  if ( order[0] < 0 ) {
218  PZError << "Compelbas::shape --> Invalid dimension for arguments: order = " << order[0]
219  << " phi.Rows = " << (int) phi.Rows() << " dphi.Cols = " << (int) dphi.Cols() << "\n";
220  return;
221  }
222  if(phi.Rows() < order[0]+1) {
223  PZError << "TPZShapeLinear::shape --> Invalid dimension for argument phi " << endl;
224  phi.Resize(order[0]+1, phi.Cols());
225  }
226  if(dphi.Cols() < order[0]+1) {
227  PZError << "TPZShapeLinear::shape --> Invalid dimension for argument dphi " << endl;
228  dphi.Resize(dphi.Rows(),order[0]+1);
229  }
230 #endif
231 
232  if ( order[0] == 0)
233  {
234  phi(0,0) = 1.;
235  dphi(0,0) = 0.;
236  } else
237  { // Linear shape functions
238  phi(0,0) = (1-x[0])/2.;
239  phi(1,0) = (1+x[0])/2.;
240  dphi(0,0) = -0.5;
241  dphi(0,1)= 0.5;
242  }
243 
244  int is,d;
245  TPZFNMatrix<100> phiblend(NSides,1),dphiblend(Dimension,NSides);
246  for(is=0; is<NCornerNodes; is++)
247  {
248  phiblend(is,0) = phi(is,0);
249  for(d=0; d<Dimension; d++)
250  {
251  dphiblend(d,is) = dphi(d,is);
252  }
253  }
254  ShapeGenerating(x,phiblend,dphiblend);
255  // Quadratic or higher shape functions
256  int num2 = order[0]-1;
257  int transformationindex = GetTransformId1d(id);
258  TPZFNMatrix<10> phiint(num2,1),dphiint(1,num2);
259  if(num2 > 0)
260  {
261  ShapeInternal(x,order[0],phiint,dphiint,transformationindex);
262  }
263  int ord;
264  for (ord = 2; ord < order[0]+1; ord++) { // even functions
265  dphi(0,ord) = dphiint(0,ord-2)*phiblend(2,0)+dphiblend(0,2)*phiint(ord-2,0);
266  phi(ord,0) = phiint(ord-2,0)*phiblend(2,0);
267  }
268  }
269  void TPZShapeLinear::ShapeCorner(TPZVec<REAL> &pt,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi){
270 
271  }
272  void TPZShapeLinear::SideShape(int side, TPZVec<REAL> &pt, TPZVec<int64_t> &id, TPZVec<int> &order,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
273  switch(side) {
274  case 0:
275  case 1:
276  TPZShapePoint::Shape(pt,id,order,phi,dphi);
277  break;
278  case 2:
279  Shape(pt,id,order,phi,dphi);
280  break;
281  }
282  }
283 
284  void TPZShapeLinear::ShapeOrder(TPZVec<int64_t> &id, TPZVec<int> &order, TPZGenMatrix<int> &shapeorders)//, TPZVec<int64_t> &sides
285  {
286  int nshape = 2+(order[0]-1);
287  if (shapeorders.Rows() != nshape) {
288  DebugStop();
289  }
290  shapeorders(0,0) = 1;
291  shapeorders(1,0) = 1;
292  for (int i=2; i<nshape; i++) {
293  shapeorders(i,0) = i;
294  }
295  }
296 
297 
298  void TPZShapeLinear::SideShapeOrder(int side, TPZVec<int64_t> &id, int order, TPZGenMatrix<int> &shapeorders)
299  {
300  DebugStop();
301  }
302  //versao nova-> o ponto vem transformado
303  void TPZShapeLinear::ShapeInternal(TPZVec<REAL> &x, int ord,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi){
304  // Quadratic or higher shape functions
305  int num = ord-1;
306  if(num <= 0) return;
307  REAL y;
308  fOrthogonal(x[0],num,phi,dphi);
309  }
310 
311  //versao antiga-> o ponto precisa ser transformado
312  void TPZShapeLinear::ShapeInternal(TPZVec<REAL> &x, int ord,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,int transformation_index){
313  // Quadratic or higher shape functions
314  int num = ord-1;
315  if(num <= 0) return;
316  REAL y;
317  TransformPoint1d(transformation_index, x[0], y);
318  fOrthogonal(y,num,phi,dphi);
319  TransformDerivative1d(transformation_index, num, dphi);
320  }
321 
322 
323 
324  void TPZShapeLinear::ShapeInternal(int side, TPZVec<REAL> &x, int order,TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi){
325  if (side != 2) {
326  return;
327  }
328  ShapeInternal(x, order, phi, dphi);
329  }
330 
331  void TPZShapeLinear::TransformPoint1d(int transid,REAL in,REAL &out) {
332  if (!transid) out = in;
333  else out = -in;
334  }
335  TPZTransform<REAL> TPZShapeLinear::ParametricTransform(int transid) {
336  TPZTransform<REAL> trans(1,1);
337  if (!transid) trans.Mult()(0,0) = 1;
338  else trans.Mult()(0,0) = -1;
339  return trans;
340  }
341 
342 
343  void TPZShapeLinear::TransformPoint1d(int transid,double &val) {
344 
345  if (!transid) val = 1.0;
346  else val = -1.0;
347  }
348 
349 
350  void TPZShapeLinear::TransformDerivative1d(int transid,int num,TPZFMatrix<REAL> &in) {
351 
352  if(transid == 0) return;
353  int i;
354  for(i=0;i<num;i++) in(0,i) = -in(0,i);
355  }
356 
357  int TPZShapeLinear::GetTransformId1d(TPZVec<int64_t> &id) {
358  if (id[1] < id[0]) return 1;
359  else return 0;
360  }
361 
362  int TPZShapeLinear::NConnectShapeF(int side, int order) {
363  if(side<2) return 1;//0 a 4
364  if(side<3) return (order-1);//6 a 14
365  PZError << "TPZShapeLinear::NConnectShapeF, bad parameter side " << side << endl;
366  return 0;
367  }
368 
369  int TPZShapeLinear::NShapeF(TPZVec<int> &order) {
370  int in,res=NCornerNodes;
371  for(in=NCornerNodes;in<NSides;in++) res += NConnectShapeF(in,order[in-NCornerNodes]);
372  return res;
373  }
374 
375 #ifdef _AUTODIFF
376  void TPZShapeLinear::ShapeInternal(FADREAL & x,int num,TPZVec<FADREAL> & phi,int transformation_index){
377  // Quadratic or higher shape functions
378  if(num <= 0) return;
379  FADREAL y;
380  TransformPoint1d(transformation_index,x,y);
381  FADfOrthogonal(y,num,phi);
382  // TransformDerivative1d(transformation_index,num,phi);
383  }
384 
385  void TPZShapeLinear::TransformPoint1d(int transid,FADREAL & in,FADREAL &out) {
386  if (!transid) out = in;
387  else out = -in;
388  }
389 
390  void TPZShapeLinear::Chebyshev(FADREAL & x,int num,TPZVec<FADREAL> &phi){
391  // Quadratic or higher shape functions
392  if(num <= 0) return;
393  //phi.Put(0,0,1.0);
394  //dphi.Put(0,0, 0.0);
395  phi[0] = 1.0; // <!> Remark: the derivatives other than the 0th are set to null
396  if(num == 1) return;
397  //phi.Put(1,0, x);
398  //dphi.Put(0,1, 1.0);
399  phi[1] = x;
400  //phi[1].fastAccessDx(0)=1.0; // <!> Remark: the derivatives other than the 0th aren't set to null
401  //just ensuring the derivatives are properly initialized and that FAD objects of more than
402  // one derivative are used
403  int ord;
404  for(ord = 2;ord<num;ord++) {
405  //phi.Put(ord,0, 2.0*x*phi(ord-1,0) - phi(ord-2,0));
406  //dphi.Put(0,ord, 2.0*x*dphi(0,ord-1) + 2.0*phi(ord-1,0) - dphi(0,ord-2));
407  phi[ord] = x * phi[ord-1] * 2.0 - phi[ord-2];
408  }
409  }
410 
411  void (*TPZShapeLinear::FADfOrthogonal)(FADREAL&,int ,TPZVec<FADREAL> &) = TPZShapeLinear::Chebyshev/*(FADREAL&, int, TPZVec<FADREAL>&)*/;//Chebyshev;
412 
413 #endif
414 
415 };
pzshape
groups all classes dedicated to the computation of shape functions
Definition: pzshapeextend.h:16
TPZTransform::Mult
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
pzerror.h
Defines PZError.
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
pzshapelinear.h
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
Shape
void Shape(TPZVec< REAL > &pt, TPZVec< int > orders, TPZVec< TPZTransform< REAL > > &transvec, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi)
Definition: pzgenericshape.h:32
TPZMatrix::Put
virtual int Put(const int64_t row, const int64_t col, const TVar &value)
Put values with bounds checking if DEBUG variable is defined.
Definition: pzmatrix.h:836
test.res
string res
Definition: test.py:151
pzshapepoint.h
Contains TPZShapePoint class which implements the shape function associated with a point...
TPZGenMatrix
Implements generic class which holds a matrix of objects. Matrix.
Definition: pzshtmat.h:18
TPZGenMatrix::Rows
int64_t Rows() const
Definition: pzshtmat.h:45
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
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
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
Generated on Sun Aug 16 2020 11:17:09 for NeoPZ by   doxygen 1.8.13