NeoPZ
convtest.cpp
Go to the documentation of this file.
1 
5 #include "convtest.h"
6 
7 ConvTest::ConvTest()
8 {
9 }
10 
11 
12 ConvTest::~ConvTest()
13 {
14 }
15 
16 #include "pzgeoel.h"
17 
19 void ConvTest::JacobianConv(TPZGeoEl &Object, TPZVec< REAL > StartPoint)
20 {
21  TPZVec< REAL > Out(3,0.);
22  if(Object.Dimension() == 1)
23  {
24  TPZFMatrix<REAL> jacobian(1,1);
25  TPZFMatrix<REAL> Axes(1,3);
26  REAL detJacobian;
27  TPZFMatrix<REAL> InvJac(1,1);
28  TPZVec< REAL > QsiEtaIni(1);
29  QsiEtaIni[0] = StartPoint[0];
30  const double deltaQsi = 0.1;
31  double alpha;
32 
33  std::cout << "\ninitial Qsi = " << QsiEtaIni[0] << "\n";
34  std::cout << "deltaQsi = const = " << deltaQsi << "\n\n";
35 
36  TPZVec< REAL > OutAprox(3);
37  TPZFMatrix<REAL> error(11,1,0.);
38  double dX, dY, dZ;
39 
40  Object.Jacobian(QsiEtaIni,jacobian,Axes,detJacobian,InvJac);
41 
42  for(int i = 0; i < 11; i++)
43  {
44  alpha = i/10.;
45  Object.X(QsiEtaIni,Out);//for aproximate compute
46  // std::cout << "\nalpha = " << alpha << std::endl;
47  // std::cout << "f(x) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
48 
49  dX = alpha*( jacobian.GetVal(0,0)*deltaQsi )*Axes(0,0);
50  OutAprox[0] = Out[0] + dX;
51 
52  dY = alpha*( jacobian.GetVal(0,0)*deltaQsi )*Axes(0,1);
53  OutAprox[1] = Out[1] + dY;
54 
55  dZ = alpha*( jacobian.GetVal(0,0)*deltaQsi )*Axes(0,2);
56  OutAprox[2] = Out[2] + dZ;
57 
58  StartPoint[0] = QsiEtaIni[0] + alpha*deltaQsi;
59 
60  Object.X(StartPoint,Out);//for real compute
61 
62  // std::cout << "alpha.(axes.J).dx : x = " << dX << " | y = " << dY << " | z = " << dZ << "\n";
63  // std::cout << "f(x) + alpha.(axes.J).dx : x = " << OutAprox[0] << " | y = " << OutAprox[1] << " | z = " << OutAprox[2] << "\n";
64  // std::cout << "f(x + alpha.dx) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
65  // std::cout << "--------------------------------------------------------------------\n";
66 
67  Out[0] -= OutAprox[0];
68  Out[1] -= OutAprox[1];
69  Out[2] -= OutAprox[2];
70 
71  double XDiffNorm = sqrt(Out[0]*Out[0] + Out[1]*Out[1] + Out[2]*Out[2]);
72  error(int(i),0) = XDiffNorm;
73  }
74 
75  // std::cout << "ERROR Vector:\n"; error.Print();
76 
77  std::cout << "Convergence Order:\n";
78  for(int j = 2; j < 11; j++)
79  {
80  std::cout << ( log(error(1,0)) - log(error(j,0)) )/( log(0.1)-log(j/10.) ) << "\n";
81  }
82  std::cout << "\nIf another kind of results are needed, edit the ConvTest Class on source code!\n";
83  return;
84  }
85 
86  if(Object.Dimension() == 2)
87  {
88  TPZFMatrix<REAL> jacobian(2,2);
89  TPZFMatrix<REAL> Axes(2,3);
90  REAL detJacobian;
91  TPZFMatrix<REAL> InvJac(2,2);
92  TPZVec< REAL > QsiEtaIni(2);
93  QsiEtaIni[0] = StartPoint[0];
94  QsiEtaIni[1] = StartPoint[1];
95  const double deltaQsi = 0.1;
96  const double deltaEta = 0.1;
97  double alpha;
98 
99  std::cout << "\ninitial Qsi = " << QsiEtaIni[0] << " | initial Eta = " << QsiEtaIni[1] << "\n";
100  std::cout << "deltaQsi = const = " << deltaQsi << " | deltaEta = const = " << deltaEta << "\n\n";
101 
102  TPZVec< REAL > OutAprox(3);
103  TPZFMatrix<REAL> error(11,1,0.);
104  double dX, dY, dZ;
105 
106  Object.Jacobian(QsiEtaIni,jacobian,Axes,detJacobian,InvJac);
107 
108  for(int i = 0; i <= 10; i++)
109  {
110  alpha = i/10.;
111  Object.X(QsiEtaIni,Out);//for aproximate compute
112  // std::cout << "\nalpha = " << alpha << std::endl;
113  // std::cout << "f(x) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
114 
115  dX = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta )*Axes(0,0) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta )*Axes(1,0);
116  OutAprox[0] = Out[0] + dX;
117 
118  dY = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta )*Axes(0,1) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta )*Axes(1,1);
119  OutAprox[1] = Out[1] + dY;
120 
121  dZ = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta )*Axes(0,2) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta )*Axes(1,2);
122  OutAprox[2] = Out[2] + dZ;
123 
124  StartPoint[0] = QsiEtaIni[0] + alpha*deltaQsi;
125  StartPoint[1] = QsiEtaIni[1] + alpha*deltaEta;
126  Object.X(StartPoint,Out);//for real compute
127 
128  // std::cout << "alpha.(axes.J).dx : x = " << dX << " | y = " << dY << " | z = " << dZ << "\n";
129  // std::cout << "f(x) + alpha.(axes.J).dx : x = " << OutAprox[0] << " | y = " << OutAprox[1] << " | z = " << OutAprox[2] << "\n";
130  // std::cout << "f(x + alpha.dx) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
131  // std::cout << "--------------------------------------------------------------------\n";
132 
133  Out[0] -= OutAprox[0];
134  Out[1] -= OutAprox[1];
135  Out[2] -= OutAprox[2];
136 
137  double XDiffNorm = sqrt(Out[0]*Out[0] + Out[1]*Out[1] + Out[2]*Out[2]);
138  error(int(i),0) = XDiffNorm;
139  }
140 
141  // std::cout << "ERROR Vector:\n"; error.Print();
142 
143  std::cout << "Convergence Order:\n";
144  for(int j = 2; j < 11; j++)
145  {
146  std::cout << ( log(error(1,0)) - log(error(j,0)) )/( log(0.1)-log(j/10.) ) << "\n";
147  }
148  std::cout << "\nIf another kind of results are needed, edit the ConvTest Class on source code!\n";
149  return;
150  }
151 
152  if(Object.Dimension() == 3)
153  {
154  TPZFMatrix<REAL> jacobian(3,3);
155  TPZFMatrix<REAL> Axes(3,3);
156  REAL detJacobian;
157  TPZFMatrix<REAL> InvJac(3,3);
158  TPZVec< REAL > QsiEtaIni(3);
159  QsiEtaIni[0] = StartPoint[0];
160  QsiEtaIni[1] = StartPoint[1];
161  QsiEtaIni[2] = StartPoint[2];
162  const double deltaQsi = 0.1;
163  const double deltaEta = 0.1;
164  const double deltaZeta = 0.1;
165  double alpha;
166 
167  std::cout << "\ninitial Qsi = " << QsiEtaIni[0] << " | initial Eta = " << QsiEtaIni[1] << " | initial Zeta = " << QsiEtaIni[2] <<"\n";
168  std::cout << "deltaQsi = const = " << deltaQsi << " | deltaEta = const = " << deltaEta << " | deltaZeta = const = " << deltaZeta <<"\n\n";
169 
170  TPZVec< REAL > OutAprox(3);
171  TPZFMatrix<REAL> error(11,1,0.);
172  double dX, dY, dZ;
173 
174  Object.Jacobian(QsiEtaIni,jacobian,Axes,detJacobian,InvJac);
175 
176  for(int i = 0; i <= 10; i++)
177  {
178  alpha = i/10.;
179  Object.X(QsiEtaIni,Out);//for aproximate compute
180  // std::cout << "\nalpha = " << alpha << std::endl;
181  // std::cout << "f(x) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
182 
183  dX = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta + jacobian.GetVal(0,2)*deltaZeta)*Axes(0,0) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta + jacobian.GetVal(1,2)*deltaZeta)*Axes(1,0) + alpha*( jacobian.GetVal(2,0)*deltaQsi + jacobian.GetVal(2,1)*deltaEta + jacobian.GetVal(2,2)*deltaZeta)*Axes(2,0);
184  OutAprox[0] = Out[0] + dX;
185 
186  dY = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta + jacobian.GetVal(0,2)*deltaZeta)*Axes(0,1) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta + jacobian.GetVal(1,2)*deltaZeta)*Axes(1,1) + alpha*( jacobian.GetVal(2,0)*deltaQsi + jacobian.GetVal(2,1)*deltaEta + jacobian.GetVal(2,2)*deltaZeta)*Axes(2,1);
187  OutAprox[1] = Out[1] + dY;
188 
189  dZ = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta + jacobian.GetVal(0,2)*deltaZeta)*Axes(0,2) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta + jacobian.GetVal(1,2)*deltaZeta)*Axes(1,2) + alpha*( jacobian.GetVal(2,0)*deltaQsi + jacobian.GetVal(2,1)*deltaEta + jacobian.GetVal(2,2)*deltaZeta)*Axes(2,2);
190  OutAprox[2] = Out[2] + dZ;
191 
192  StartPoint[0] = QsiEtaIni[0] + alpha*deltaQsi;
193  StartPoint[1] = QsiEtaIni[1] + alpha*deltaEta;
194  StartPoint[2] = QsiEtaIni[2] + alpha*deltaZeta;
195  Object.X(StartPoint,Out);//for real compute
196 
197  // std::cout << "alpha.(axes.J).dx : x = " << dX << " | y = " << dY << " | z = " << dZ << "\n";
198  // std::cout << "f(x) + alpha.(axes.J).dx : x = " << OutAprox[0] << " | y = " << OutAprox[1] << " | z = " << OutAprox[2] << "\n";
199  // std::cout << "f(x + alpha.dx) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
200  // std::cout << "--------------------------------------------------------------------\n";
201 
202  Out[0] -= OutAprox[0];
203  Out[1] -= OutAprox[1];
204  Out[2] -= OutAprox[2];
205 
206  double XDiffNorm = sqrt(Out[0]*Out[0] + Out[1]*Out[1] + Out[2]*Out[2]);
207  error(int(i),0) = XDiffNorm;
208  }
209 
210  // std::cout << "ERROR Vector:\n"; error.Print();
211 
212  std::cout << "Convergence Order:\n";
213  for(int j = 2; j < 11; j++)
214  {
215  std::cout << ( log(error(1,0)) - log(error(j,0)) )/( log(0.1)-log(j/10.) ) << "\n";
216  }
217  std::cout << "\nIf another kind of results are needed, edit the ConvTest Class on source code!\n";
218  return;
219  }
220  else
221  {
222  std::cout << "Element don't fit in an option of Convergence Analys!\nSee ConvTest Class...\n";
223  return;
224  }
225 }
226 
228 void ConvTest::JacobianConv(TPZGeoElSide &Object, TPZVec< REAL > StartPoint)
229 {
230  TPZVec< REAL > Out(3,0.);
231  if(Object.Dimension() == 1)
232  {
233  TPZFMatrix<REAL> jacobian(1,1);
234  TPZFMatrix<REAL> Axes(1,3);
235  REAL detJacobian;
236  TPZFMatrix<REAL> InvJac(1,1);
237  TPZVec< REAL > QsiEtaIni(1);
238  QsiEtaIni[0] = StartPoint[0];
239  const double deltaQsi = 0.1;
240  double alpha;
241 
242  std::cout << "\ninitial Qsi = " << QsiEtaIni[0] << "\n";
243  std::cout << "deltaQsi = const = " << deltaQsi << "\n\n";
244 
245  TPZVec< REAL > OutAprox(3);
246  TPZFMatrix<REAL> error(11,1,0.);
247  double dX, dY, dZ;
248 
249  Object.Jacobian(QsiEtaIni,jacobian,Axes,detJacobian,InvJac);
250 
251  for(int i = 0; i < 11; i++)
252  {
253  alpha = i/10.;
254  Object.X(QsiEtaIni,Out);//for aproximate compute
255  // std::cout << "\nalpha = " << alpha << std::endl;
256  // std::cout << "f(x) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
257 
258  dX = alpha*( jacobian.GetVal(0,0)*deltaQsi )*Axes(0,0);
259  OutAprox[0] = Out[0] + dX;
260 
261  dY = alpha*( jacobian.GetVal(0,0)*deltaQsi )*Axes(0,1);
262  OutAprox[1] = Out[1] + dY;
263 
264  dZ = alpha*( jacobian.GetVal(0,0)*deltaQsi )*Axes(0,2);
265  OutAprox[2] = Out[2] + dZ;
266 
267  StartPoint[0] = QsiEtaIni[0] + alpha*deltaQsi;
268 
269  Object.X(StartPoint,Out);//for real compute
270 
271  // std::cout << "alpha.(axes.J).dx : x = " << dX << " | y = " << dY << " | z = " << dZ << "\n";
272  // std::cout << "f(x) + alpha.(axes.J).dx : x = " << OutAprox[0] << " | y = " << OutAprox[1] << " | z = " << OutAprox[2] << "\n";
273  // std::cout << "f(x + alpha.dx) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
274  // std::cout << "--------------------------------------------------------------------\n";
275 
276  Out[0] -= OutAprox[0];
277  Out[1] -= OutAprox[1];
278  Out[2] -= OutAprox[2];
279 
280  double XDiffNorm = sqrt(Out[0]*Out[0] + Out[1]*Out[1] + Out[2]*Out[2]);
281  error(int(i),0) = XDiffNorm;
282  }
283 
284  // std::cout << "ERROR Vector:\n"; error.Print();
285 
286  std::cout << "Convergence Order:\n";
287  for(int j = 2; j < 11; j++)
288  {
289  std::cout << ( log(error(1,0)) - log(error(j,0)) )/( log(0.1)-log(j/10.) ) << "\n";
290  }
291  std::cout << "\nIf another kind of results are needed, edit the ConvTest Class on source code!\n";
292  return;
293  }
294 
295  if(Object.Dimension() == 2)
296  {
297  TPZFMatrix<REAL> jacobian(2,2);
298  TPZFMatrix<REAL> Axes(2,3);
299  REAL detJacobian;
300  TPZFMatrix<REAL> InvJac(2,2);
301  TPZVec< REAL > QsiEtaIni(2);
302  QsiEtaIni[0] = StartPoint[0];
303  QsiEtaIni[1] = StartPoint[1];
304  const double deltaQsi = 0.1;
305  const double deltaEta = 0.1;
306  double alpha;
307 
308  std::cout << "\ninitial Qsi = " << QsiEtaIni[0] << " | initial Eta = " << QsiEtaIni[1] << "\n";
309  std::cout << "deltaQsi = const = " << deltaQsi << " | deltaEta = const = " << deltaEta << "\n\n";
310 
311  TPZVec< REAL > OutAprox(3);
312  TPZFMatrix<REAL> error(11,1,0.);
313  double dX, dY, dZ;
314 
315  Object.Jacobian(QsiEtaIni,jacobian,Axes,detJacobian,InvJac);
316 
317  for(int i = 0; i <= 10; i++)
318  {
319  alpha = i/10.;
320  Object.X(QsiEtaIni,Out);//for aproximate compute
321  // std::cout << "\nalpha = " << alpha << std::endl;
322  // std::cout << "f(x) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
323 
324  dX = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta )*Axes(0,0) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta )*Axes(1,0);
325  OutAprox[0] = Out[0] + dX;
326 
327  dY = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta )*Axes(0,1) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta )*Axes(1,1);
328  OutAprox[1] = Out[1] + dY;
329 
330  dZ = alpha*( jacobian.GetVal(0,0)*deltaQsi + jacobian.GetVal(0,1)*deltaEta )*Axes(0,2) + alpha*( jacobian.GetVal(1,0)*deltaQsi + jacobian.GetVal(1,1)*deltaEta )*Axes(1,2);
331  OutAprox[2] = Out[2] + dZ;
332 
333  StartPoint[0] = QsiEtaIni[0] + alpha*deltaQsi;
334  StartPoint[1] = QsiEtaIni[1] + alpha*deltaEta;
335  Object.X(StartPoint,Out);//for real compute
336 
337  // std::cout << "alpha.(axes.J).dx : x = " << dX << " | y = " << dY << " | z = " << dZ << "\n";
338  // std::cout << "f(x) + alpha.(axes.J).dx : x = " << OutAprox[0] << " | y = " << OutAprox[1] << " | z = " << OutAprox[2] << "\n";
339  // std::cout << "f(x + alpha.dx) : x = " << Out[0] << " | y = " << Out[1] << " | z = " << Out[2] << "\n";
340  // std::cout << "--------------------------------------------------------------------\n";
341 
342  Out[0] -= OutAprox[0];
343  Out[1] -= OutAprox[1];
344  Out[2] -= OutAprox[2];
345 
346  double XDiffNorm = sqrt(Out[0]*Out[0] + Out[1]*Out[1] + Out[2]*Out[2]);
347  error(int(i),0) = XDiffNorm;
348  }
349 
350  // std::cout << "ERROR Vector:\n"; error.Print();
351 
352  std::cout << "Convergence Order:\n";
353  for(int j = 2; j < 11; j++)
354  {
355  std::cout << ( log(error(1,0)) - log(error(j,0)) )/( log(0.1)-log(j/10.) ) << "\n";
356  }
357  std::cout << "\nIf another kind of results are needed, edit the ConvTest Class on source code!\n";
358  return;
359  }
360 
361  else
362  {
363  std::cout << "Element don't fit in an option of Convergence Analys!\nSee ConvTest Class...\n";
364  return;
365  }
366 }
ConvTest::JacobianConv
void JacobianConv(TPZGeoEl &Object, TPZVec< REAL > QsiEta)
Evaluates the Jacobian by obtained Convergence Order.
Definition: convtest.cpp:19
ConvTest::ConvTest
ConvTest()
Default constructor.
Definition: convtest.cpp:7
convtest.h
Contains ConvTest class which implements methods to evaluate jacobians by obtained convergence order ...
TPZGeoElSide
Utility class which represents an element with its side. The Geometric approximation classes Geometry...
Definition: pzgeoelside.h:83
error
void error(char *string)
Definition: testShape.cc:7
TPZGeoEl
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
TPZGeoElSide::X
void X(TPZVec< REAL > &loc, TPZVec< REAL > &result) const
X coordinate of a point loc of the side.
Definition: pzgeoelside.cpp:117
TPZGeoEl::Jacobian
void Jacobian(TPZVec< REAL > &qsi, TPZFMatrix< REAL > &jac, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Compute a decomposition of the gradient of the mapping function, as a rotation matrix (Jacobian) and ...
Definition: pzgeoel.cpp:1144
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
TPZGeoEl::Dimension
virtual int Dimension() const =0
Returns the dimension of the element.
log
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_ log
Definition: tfadfunc.h:130
TPZGeoEl::X
virtual void X(TPZVec< REAL > &qsi, TPZVec< REAL > &result) const =0
Return the coordinate in real space of the point coordinate in the master element space...
TPZGeoElSide::Dimension
int Dimension() const
the dimension associated with the element/side
Definition: pzgeoelside.cpp:593
TPZGeoElSide::Jacobian
void Jacobian(TPZVec< REAL > &param, TPZFMatrix< REAL > &jacobian, TPZFMatrix< REAL > &axes, REAL &detjac, TPZFMatrix< REAL > &jacinv) const
Jacobian associated with the side of the element.
Definition: pzgeoelside.cpp:189
ConvTest::~ConvTest
~ConvTest()
Default destructor.
Definition: convtest.cpp:12
TPZFMatrix::GetVal
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
Generated on Sun Aug 16 2020 11:17:10 for NeoPZ by   doxygen 1.8.13