NeoPZ
pzshapetriang.cpp
Go to the documentation of this file.
1 
6 #include "pzshapetriang.h"
7 #include "pzshapelinear.h"
8 #include "pzshapepoint.h"
9 #include "pzmanvector.h"
10 #include "pzerror.h"
11 #include "pzreal.h"
12 #include "pzgenericshape.h"
13 
14 using namespace std;
15 
16 namespace pzshape {
17 
19  REAL TPZShapeTriang::gTrans2dT[6][2][2] = {//s* , t*
20  { { 1., 0.},{ 0., 1.} },
21  { { 0., 1.},{ 1., 0.} },
22  { { 0., 1.},{-1.,-1.} },//s* = t t* = -s-t-1 , etc
23  { {-1.,-1.},{ 0., 1.} },
24  { {-1.,-1.},{ 1., 0.} },
25  { { 1., 0.},{-1.,-1.} }
26  };
27 
28  REAL TPZShapeTriang::gVet2dT[6][2] = { {0.,0.},{0.,0.},{0.,1.},{1.,0.},{1.,0.},{0.,1.} };
29 
30  REAL TPZShapeTriang::gRibTrans2dT1d[3][2] = { {2.,1.},{-1.,1.},{-1.,-2.} };//Cedric : 06/03/99
31 
32  REAL TPZShapeTriang::gVet1dT[3] = {-1.,0.,1.};
33 
34 
35  void TPZShapeTriang::ShapeCorner(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
36 
37  phi(0,0) = 1.-pt[0]-pt[1];
38  phi(1,0) = pt[0];
39  phi(2,0) = pt[1];
40  dphi(0,0) = -1.;
41  dphi(1,0) = -1.;
42  dphi(0,1) = 1.;
43  dphi(1,1) = 0.;
44  dphi(0,2) = 0.;
45  dphi(1,2) = 1.;
46  }
47 
54  void TPZShapeTriang::ShapeGenerating(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi)
55  {
56  int is;
57  for(is=3; is<6; is++)
58  {
59  int is1 = is%3;
60  int is2 = (is+1)%3;
61  phi(is,0) = phi(is1,0)*phi(is2,0);
62  dphi(0,is) = dphi(0,is1)*phi(is2,0)+phi(is1,0)*dphi(0,is2);
63  dphi(1,is) = dphi(1,is1)*phi(is2,0)+phi(is1,0)*dphi(1,is2);
64  }
65  int is1 = 0;
66  int is2 = 1;
67  int is3 = 2;
68  phi(is,0) = phi(is1,0)*phi(is2,0)*phi(is3,0);
69  dphi(0,is) = dphi(0,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(0,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(0,is3);
70  dphi(1,is) = dphi(1,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(1,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(1,is3);
71 
72  // Make the generating shape functions linear and unitary
73  REAL mult[] = {1.,1.,1.,4.,4.,4.,27.};
74  for(is=3;is<NSides; is++)
75  {
76  phi(is,0) *= mult[is];
77  dphi(0,is) *= mult[is];
78  dphi(1,is) *= mult[is];
79  }
80 
81 
82  }
83 
84  void TPZShapeTriang::Shape(TPZVec<REAL> &pt, TPZVec<int64_t> &id, TPZVec<int> &order,
85  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
86  ShapeCorner(pt,phi,dphi);
87  if (order[0] < 2 && order[1] < 2 && order[2] < 2 && order[3] < 3) return;
88  int is,d;
89 
90  TPZFNMatrix<50,REAL> phiblend(NSides,1);
91  TPZFNMatrix<100,REAL> dphiblend(Dimension,NSides);
92  for(is=0; is<NCornerNodes; is++)
93  {
94  phiblend(is,0) = phi(is,0);
95  for(d=0; d<Dimension; d++)
96  {
97  dphiblend(d,is) = dphi(d,is);
98  }
99  }
100  ShapeGenerating(pt,phiblend,dphiblend);
101 
102  REAL out;
103  int shape = 3;
104  for (int rib = 0; rib < 3; rib++) {
105 
106  ProjectPoint2dTriangToRib(rib,pt,out);
107  TPZManVector<REAL,1> outvec(1,out);
108  TPZVec<int64_t> ids(2);
109  ids[0] = id[rib%3];
110  ids[1] = id[(rib+1)%3];
111  int ord2 = order[rib]-1;//numero de shapes por lado rib
112  TPZFNMatrix<50,REAL> phin(ord2,1);
113  TPZFNMatrix<100,REAL>dphin(2,ord2);
114  TPZShapeLinear *shplin=0;
115  shplin->ShapeInternal(outvec,order[rib],phin,dphin,shplin->GetTransformId1d(ids));
116 // dphin.Print("dphin= ",std::cout,EMathematicaInput);
117  //TransformDerivativeFromRibToTriang(rib,ord2,dphin);
118  for (int i = 0; i < ord2; i++) {
119  phi(shape,0) = phiblend(rib+3,0)*phin(i,0);
120  for(int xj=0;xj<2;xj++) {
121  dphi(xj,shape) = dphiblend(xj,rib+3)* phin( i, 0)+
122  phiblend(rib+3, 0 )* dphin(xj,i);
123  }
124  shape++;
125  }
126  }
127 
128  if (order[3] < 3) return;//ordem na face
129  int ord = order[3]-2;//num de shapes da face
130  int nsh = (ord*(ord+1))/2;
131  TPZFNMatrix<50,REAL> phin(nsh,1);
132  TPZFNMatrix<100,REAL> dphin(2,nsh);
133  ShapeInternal(pt,order[3]-2, phin, dphin,GetTransformId2dT(id));
134 // dphin.Print("dphin= ",std::cout,EMathematicaInput);
135 
136  for(int i=0;i<nsh;i++) {//number of internal shape equal maximal order
137  phi(shape,0) = phiblend(6,0)*phin(i,0);
138  for(int d=0;d<2;d++) {
139  dphi(d,shape) = dphiblend(d,6)* phin(i,0) +
140  phiblend(6,0)*dphin(d,i);
141  }
142  shape++;
143  }
144  }
145 
146  void TPZShapeTriang::SideShape(int side, TPZVec<REAL> &pt, TPZVec<int64_t> &id, TPZVec<int> &order,
147  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
148  if(side<0 || side>6) PZError << "TPZShapeTriang::SideShape. Bad paramenter side.\n";
149  else if(side==6) Shape(pt,id,order,phi,dphi);
150  else if(side<3) {
151  TPZShapePoint::Shape(pt,id,order,phi,dphi);
152  } else {
153  TPZShapeLinear::Shape(pt,id,order,phi,dphi);
154  }
155 
156 
157  }
158 
159  void TPZShapeTriang::ShapeOrder(TPZVec<int64_t> &id, TPZVec<int> &order, TPZGenMatrix<int> &shapeorders)//, TPZVec<int64_t> &sides
160  {
161  int64_t nsides = TPZTriangle::NSides;
162  // o que eh o vetor order?
163  // Eu suponho que em cada posicao tem a ordem de cada lado.
164  // Na shape ja esta associado a lados com dimensao maior que 1, order[0](lado 3) ...
165 
166  int nshape = NCornerNodes;
167 
168  for (int is = NCornerNodes; is < nsides; is++)
169  {
170  nshape += NConnectShapeF(is,order[is-NCornerNodes]);
171  }
172  // shapeorders ja vem com tamanho definido, entao a conta acima so serve para ver se as ordens estao corretas
173  if (nshape!=shapeorders.Rows())
174  {
175  // Algo errado
176  DebugStop();
177  }
178 
179  int linha = 0;
180  for (int side = 0; side < nsides; side++)
181  {
182 
183  //nshape = NConnectShapeF(side,order[side]);
184  //int sideorder = order[side];
185  nshape = 1;
186  if(side >= NCornerNodes) nshape = NConnectShapeF(side,order[side-NCornerNodes]);
187  int sideorder = 1;
188  if(side >= NCornerNodes) sideorder = order[side-NCornerNodes];
189 
190  TPZGenMatrix<int> locshapeorders(nshape,3);
191  SideShapeOrder(side, id, sideorder, locshapeorders);
192 
193  int nlin = locshapeorders.Rows();
194  int ncol = locshapeorders.Cols();
195 
196  for (int il = 0; il<nlin; il++)
197  {
198  for (int jc = 0; jc<ncol; jc++)
199  {
200  shapeorders(linha, jc) = locshapeorders(il, jc);
201  }
202  linha++;
203  }
204  }
205 
206  }
207 
208 
209  void TPZShapeTriang::SideShapeOrder(int side, TPZVec<int64_t> &id, int order, TPZGenMatrix<int> &shapeorders)
210  {
211 
212  if (side<=2)
213  {
214  if(shapeorders.Rows() != 1)
215  {
216  DebugStop();
217  }
218  shapeorders(0,0) = 1;
219  shapeorders(0,1) = 0;
220  shapeorders(0,2) = 0;
221  }
222  else if (side == 6)
223  {
224  int nshape = (order-2)*(order-1)/2;
225  if(shapeorders.Rows() != nshape) DebugStop();
226  int nsh=1;
227  int ish = 0;
228  for (int i=3; i<=order; i++) {
229  for (int j=0; j<nsh; j++) {
230  shapeorders(ish,0) = i;
231  shapeorders(ish,1) = i;
232  ish++;
233  }
234  nsh++;
235  }
236  }
237  else
238  {
239  int nshape = order-1;
240  if (shapeorders.Rows() != nshape) {
241  DebugStop();
242  }
243  for (int ioy = 0; ioy < order-1; ioy++)
244  {
245  shapeorders(ioy,0) = ioy+2;
246  }
247  }
248 
249  }
250 
251  void TPZShapeTriang::ShapeInternal(TPZVec<REAL> &x, int order,TPZFMatrix<REAL> &phi,
252  TPZFMatrix<REAL> &dphi) {
253 
254  if((order - 2 ) <= 0) return;
255  int numshape = ((order-2)*(order-1))/2;
256 
257  TPZManVector<REAL,2> out(2,0.0);
258  out[0] = 2.*x[0]-1.;
259  out[1] = 2.*x[1]-1.;
260 
261  if (phi.Rows() < numshape || dphi.Cols() < numshape) {
262  PZError << "\nTPZCompEl::Shape2dTriangleInternal phi or dphi resized\n";
263  phi.Resize(numshape,1);
264  dphi.Resize(dphi.Rows(),numshape);
265  }
266 
267  TPZFNMatrix<50,REAL> phi0(numshape,1),phi1(numshape,1);
268  TPZFNMatrix<100,REAL> dphi0(1,numshape),dphi1(1,numshape);
269 
270  TPZShapeLinear::fOrthogonal(out[0],numshape,phi0,dphi0);
271  TPZShapeLinear::fOrthogonal(out[1],numshape,phi1,dphi1);
272  int index = 0;
273  int i;
274 
275  for (int iplusj=0;iplusj<(order - 2);iplusj++) {
276  for (int j=0;j<=iplusj;j++) {
277  i = iplusj-j;
278  phi(index,0) = phi0(i,0)*phi1(j,0);
279  dphi(0,index) = 2.0*dphi0(0,i)*phi1(j,0);
280  dphi(1,index) = 2.0*phi0(i,0)*dphi1(0,j);
281  index++;
282  }
283  }
284  }
285 
286  void TPZShapeTriang::ShapeInternal(TPZVec<REAL> &x, int order,TPZFMatrix<REAL> &phi,
287  TPZFMatrix<REAL> &dphi,int triangle_transformation_index) {
288 
289  if(order < 0) return;
290  int ord1 = order;
291  int numshape = (ord1*(ord1+1))/2;
292  TPZManVector<REAL,2> out(2);
293  TransformPoint2dT(triangle_transformation_index,x,out);
294 
295  out[0] = 2.*out[0]-1.;
296  out[1] = 2.*out[1]-1.;
297 
298  if (phi.Rows() < numshape || dphi.Cols() < numshape) {
299  PZError << "\nTPZCompEl::Shape2dTriangleInternal phi or dphi resized\n";
300  phi.Resize(numshape,1);
301  dphi.Resize(dphi.Rows(),numshape);
302  }
303 
304  TPZFNMatrix<50,REAL> phi0(numshape,1),phi1(numshape,1);
305  TPZFNMatrix<100,REAL> dphi0(1,numshape),dphi1(1,numshape);
306 
307  TPZShapeLinear::fOrthogonal(out[0],numshape,phi0,dphi0);
308  TPZShapeLinear::fOrthogonal(out[1],numshape,phi1,dphi1);
309  int index = 0;
310  int i;
311  for (int iplusj=0;iplusj<ord1;iplusj++) {
312  for (int j=0;j<=iplusj;j++) {
313  i = iplusj-j;
314  phi(index,0) = phi0(i,0)*phi1(j,0);
315  dphi(0,index) = 2.*dphi0(0,i)*phi1(j,0);
316  dphi(1,index) = 2.*phi0(i,0)*dphi1(0,j);
317  index++;
318  }
319  }
320  TransformDerivative2dT(triangle_transformation_index,numshape,dphi);
321  }
322 
323 
324 
325  void TPZShapeTriang::ShapeInternal(int side, TPZVec<REAL> &x, int order,
326  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
327  if (side < 3 || side > 6) {
328  DebugStop();
329  }
330 
331  switch (side) {
332 
333  case 3:
334  case 4:
335  case 5:
336  {
337  pzshape::TPZShapeLinear::ShapeInternal(x, order, phi, dphi);
338  }
339  break;
340  case 6:
341  {
342 
343  ShapeInternal(x, order, phi, dphi);
344  }
345  break;
346  default:
347  std::cout << __PRETTY_FUNCTION__ << " Wrong side parameter side " << side << std::endl;
348  DebugStop();
349  break;
350  }
351 
352 
353 
354  }
355 
356  void TPZShapeTriang::ProjectPoint2dTriangToRib(int rib, TPZVec<REAL> &in, REAL &out) {
357 
358  out = gRibTrans2dT1d[rib][0]*in[0]+gRibTrans2dT1d[rib][1]*in[1]+gVet1dT[rib];
359  }
360 
361  TPZTransform<REAL> TPZShapeTriang::ParametricTransform(int trans_id){
362  TPZTransform<REAL> trans(2,2);
363  trans.Mult()(0,0) = gTrans2dT[trans_id][0][0];
364  trans.Mult()(0,1) = gTrans2dT[trans_id][0][1];
365  trans.Mult()(1,0) = gTrans2dT[trans_id][1][0];
366  trans.Mult()(1,1) = gTrans2dT[trans_id][1][1];
367  trans.Sum()(0,0) =gVet2dT[trans_id][0];
368  trans.Sum()(1,0) =gVet2dT[trans_id][1];
369  return trans;
370 
371  }
372 
373  void TPZShapeTriang::TransformDerivativeFromRibToTriang(int rib,int num,TPZFMatrix<REAL> &dphi) {
374 
375  for (int j = 0;j<num;j++) {
376 
377  dphi(1,j) = gRibTrans2dT1d[rib][1]*dphi(0,j);
378  dphi(0,j) = gRibTrans2dT1d[rib][0]*dphi(0,j);
379  }
380  }
381 
382  int TPZShapeTriang::GetTransformId2dT(TPZVec<int64_t> &id) {
383 
384  int id0,id1,minid;
385  id0 = (id[0] < id[1]) ? 0 : 1;
386  minid = (id[2] < id[id0]) ? 2 : id0;
387  id0 = (minid+1)%3;
388  id1 = (minid+2)%3;
389 
390  if (id[id0] < id[id1]) {//antihorario
391 
392  if (minid == 0) return 0;
393  if (minid == 1) return 2;
394  if (minid == 2) return 4;
395 
396  } else {//horario
397 
398  if (minid == 0) return 1;
399  if (minid == 1) return 3;
400  if (minid == 2) return 5;
401  }
402  return 0;
403  }
404 
405  //transf. o ponto dentro da face triangular
406  void TPZShapeTriang::TransformPoint2dT(int transid, TPZVec<REAL> &in, TPZVec<REAL> &out) {
407 
408  out[0] = gTrans2dT[transid][0][0]*in[0]+gTrans2dT[transid][0][1]*in[1]+gVet2dT[transid][0];
409  out[1] = gTrans2dT[transid][1][0]*in[0]+gTrans2dT[transid][1][1]*in[1]+gVet2dT[transid][1];
410  }
411 
412  void TPZShapeTriang::TransformDerivative2dT(int transid, int num, TPZFMatrix<REAL> &in) {
413 
414  int i;
415  for(i=0;i<num;i++) { //ds/dcsi
416  REAL aux[2];
417  aux[0] = in(0,i);
418  aux[1] = in(1,i);
419  in(0,i) = gTrans2dT[transid][0][0]*aux[0]+gTrans2dT[transid][1][0]*aux[1];
420  in(1,i) = gTrans2dT[transid][0][1]*aux[0]+gTrans2dT[transid][1][1]*aux[1];
421  }
422  }
423 
424  int TPZShapeTriang::NConnectShapeF(int side, int order) {
425  switch(side) {
426  case 0:
427  case 1:
428  case 2:
429  return 1;
430  case 3:
431  case 4:
432  case 5:
433  return order-1;
434  case 6:
435  return (order-2) < 0 ? 0 : ((order-2)*(order-1))/2;
436  default:
437  PZError << "TPZShapeTriang::NConnectShapeF, bad parameter iconnect " << side << endl;
438  return 0;
439  }
440  }
441 
442  int TPZShapeTriang::NShapeF(TPZVec<int> &order) {
443  int in,res=NCornerNodes;
444  for(in=NCornerNodes;in<NSides;in++)
445  {
446  res += NConnectShapeF(in,order[in-NCornerNodes]);
447  }
448  return res;
449  }
450 
451 };
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
pzshape::TPZShapeLinear::ShapeInternal
static void ShapeInternal(TPZVec< REAL > &x, int ord, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, int transformation_index)
Computes the values of the orthogonal shapefunctions before multiplying them by the corner shapefunct...
Definition: pzshapelinear.cpp:312
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.
pzshapelinear.h
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
pzmanvector.h
Free store vector implementation.
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
TPZTransform::Sum
const TPZFMatrix< T > & Sum() const
Definition: pztrnsform.h:66
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::Cols
int64_t Cols() const
Definition: pzshtmat.h:47
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
pzshape::TPZShapeLinear::GetTransformId1d
static int GetTransformId1d(TPZVec< int64_t > &id)
Computes the id of the transformation which will be applied on the parameter of the element ...
Definition: pzshapelinear.cpp:357
pzshapetriang.h
Contains TPZShapeTriang class which implements the shape functions of a triangular element...
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
pzshape::TPZShapeLinear
Implements the shape functions of a linear (1D) element. Shape.
Definition: pzshapelinear.h:33
Generated on Sun Aug 16 2020 11:17:10 for NeoPZ by   doxygen 1.8.13