NeoPZ
TPZQuadSphere.h
Go to the documentation of this file.
1 //
2 // TPZQuadSphere.h
3 // PZ
4 //
5 // Created by Philippe Devloo on 3/21/14.
6 //
7 //
8 
9 #ifndef __PZ__TPZQuadSphere__
10 #define __PZ__TPZQuadSphere__
11 
12 #include <iostream>
13 #include "pzgeoquad.h"
14 
15 namespace pzgeom {
16 
17  template<class GeomQuad = pzgeom::TPZGeoQuad>
18  class TPZQuadSphere : public GeomQuad
19  {
20  private:
21  REAL fR;
22  TPZVec<REAL> fxc;
23 
24  public:
25  typedef pztopology::TPZQuadrilateral Top;
27  TPZQuadSphere(TPZVec<int64_t> &nodeindexes) : GeomQuad(nodeindexes), fR(0.), fxc(3,0.)
28  {
29  }
30 
32  TPZQuadSphere() : GeomQuad(), fR(0.), fxc(3,0.)
33  {
34  }
35 
37  TPZQuadSphere(const TPZQuadSphere &cp,
38  std::map<int64_t,int64_t> & gl2lcNdMap) : GeomQuad(cp,gl2lcNdMap), fR(0.), fxc(3,0.)
39  {
40  }
41 
43  TPZQuadSphere(const TPZQuadSphere &cp) : GeomQuad(cp), fR(cp.fR), fxc(cp.fxc)
44  {
45  }
46 
48  TPZQuadSphere(const TPZQuadSphere &cp, TPZGeoMesh &) : GeomQuad(cp), fR(cp.fR), fxc(cp.fxc)
49  {
50  }
51 
52  TPZQuadSphere &operator=(const TPZQuadSphere &cp)
53  {
54  GeomQuad::operator=(cp);
55  fR = cp.fR;
56  fxc = cp.fxc;
57  return *this;
58  }
59 
61  bool IsGeoBlendEl() const;
62 
63  static bool IsLinearMapping(int side)
64  {
65  return false;
66  }
67 
68 
70  static std::string TypeName() { return "QuadSphere";}
71 
72  void SetData(const REAL &R, const TPZVec<REAL> &xc)
73  {
74  fR = R;
75  fxc = xc;
76 #ifdef PZDEBUG
77  if (R <= 0) {
78  PZError << "R must be positive!\n";
79  DebugStop();
80  }
81  if (xc.NElements() != 3) {
82  PZError << "XCenter must have 3 coordinates!\n";
83  DebugStop();
84  }
85 #endif
86  }
87 
88  /* @brief Computes the coordinate of a point given in parameter space */
89  template<class T>
90  void X(TPZFMatrix<REAL> &coord,TPZVec<T> &loc,TPZVec<T> &result) const
91  {
92  TPZManVector<T,3> xqsi(3,0.); // will store (x,y,z) from (qsi,eta)
93  TPZManVector<T,3> xqsiLxc(3,0.); // will store (x,y,z)-xc
94  GeomQuad::X(coord,loc,xqsi); // gives the map (qsi,eta) to (x,y,z)
95 
96  T norm = 0.;
97  for (int i = 0; i < 3; i++) { // Does xqsi-xc and calculates its norm
98  xqsiLxc[i] = xqsi[i] - fxc[i];
99  norm += xqsiLxc[i] * xqsiLxc[i];
100  }
101  norm = sqrt(norm);
102 
103  for (int i = 0; i < 3; i++) {
104  result[i] = fxc[i] + xqsiLxc[i] * fR / norm;
105  }
106 
107 
108  }
109 
110  template<class T>
111  void GradX(TPZFMatrix<REAL> &coord, TPZVec<T> &param, TPZFMatrix<T> &gradx) const
112  {
113  // will first do dxdqsi and then d(phi*v)/dx, finaly d(phi*v)/dx * dxdqsi, where phi = 1/norm(xqsi - xc) and v = (xqsi - xc) * fR
114 
115  TPZManVector<T,3> xqsi(3,0.); // will store (x,y,z) from (qsi,eta)
116  TPZManVector<T,3> xqsiLxc(3,0.); // will store (x,y,z)-xc
117  GeomQuad::X(coord,param,xqsi); // gives the map (qsi,eta) to (x,y,z)
118  T norm = 0.;
119  for (int i = 0; i < 3; i++) { // Does xqsi-xc and calculates its norm
120  xqsiLxc[i] = xqsi[i] - fxc[i];
121  norm += xqsiLxc[i] * xqsiLxc[i];
122  }
123  norm = sqrt(norm);
124 
125  TPZFNMatrix<6,T> dxdqsi(3,2,0.); // But it is a (3,2) matrix. It is set (3,3) because of the later products
126  GeomQuad::GradX(coord,param,dxdqsi);
127 
128  TPZFNMatrix<3,T> gradphi(3,1,0.), v(3,1,0.); // here phi = 1/norm(xqsi - xc) and v = (xqsi - xc) * fR
129  T zero = param[0]-param[0];
130  TPZFNMatrix<9,T> gradv(3,3,zero); // here v = (xqsi - xc) * fR
131  T phi = 1./norm;
132 
133  for (int i = 0; i < 3; i++) {
134  v(i,0) = xqsiLxc[i] * fR;
135  gradv(i,i) = fR+zero;
136  gradphi(i,0) = - (1. / (norm*norm*norm) ) * xqsiLxc[i];
137  }
138 
139  TPZFNMatrix <9,T> DphivDx(3,3,0.); // will store d(phi*v)/dx
140  DphivDx = TensorProd(v,gradphi) + phi*gradv;
141 
142  DphivDx.Multiply(dxdqsi, gradx);
143  }
144 
145  /* @brief Computes the jacobian of the map between the master element and deformed element */
146 // void Jacobian(const TPZGeoEl &gel,TPZVec<REAL> &param,TPZFMatrix<REAL> &jacobian,TPZFMatrix<REAL> &axes,REAL &detjac,TPZFMatrix<REAL> &jacinv) const
147 // {
148 //
149 //
150 // // will first do dxdqsi and then d(phi*v)/dx, finaly d(phi*v)/dx * dxdqsi, where phi = 1/norm(xqsi - xc) and v = (xqsi - xc) * fR
151 //
152 // TPZManVector<REAL,3> xqsi(3,0.); // will store (x,y,z) from (qsi,eta)
153 // TPZManVector<REAL,3> xqsiLxc(3,0.); // will store (x,y,z)-xc
154 // GeomQuad::X(gel,param,xqsi); // gives the map (qsi,eta) to (x,y,z)
155 // REAL norm = 0.;
156 // for (int i = 0; i < 3; i++) { // Does xqsi-xc and calculates its norm
157 // xqsiLxc[i] = xqsi[i] - fxc[i];
158 // norm += xqsiLxc[i] * xqsiLxc[i];
159 // }
160 // norm = sqrt(norm);
161 //
162 // TPZFNMatrix<6,REAL> dxdqsi(3,2,0.); // But it is a (3,2) matrix. It is set (3,3) because of the later products
163 // DebugStop();
164 // //GeomQuad::Jacobian(gel, param, jacobian, axes, detjac, jacinv); // first calculate the derivative dxdqsi (note a lot of dummies in the parameters)
165 // TPZFMatrix<REAL> axest;
166 // axes.Transpose(&axest);
167 // axest.Multiply(jacobian, dxdqsi);
168 // jacobian.Zero();
169 //
170 //
171 // TPZFNMatrix<3,REAL> gradphi(3,1,0.), v(3,1,0.); // here phi = 1/norm(xqsi - xc) and v = (xqsi - xc) * fR
172 // TPZFNMatrix<9,REAL> gradv(3,3,0.); // here v = (xqsi - xc) * fR
173 // REAL phi = 1./norm;
174 //
175 // for (int i = 0; i < 3; i++) {
176 // v(i,0) = xqsiLxc[i] * fR;
177 // gradv(i,i) = fR;
178 // gradphi(i,0) = - (1. / (norm*norm*norm) ) * xqsiLxc[i];
179 // }
180 //
181 // TPZFNMatrix <9,REAL> DphivDx(3,3,0.); // will store d(phi*v)/dx
182 // DphivDx = TensorProd(v,gradphi) + phi*gradv;
183 //
184 // TPZFNMatrix <6,REAL> GradX(3,2,0.); // stores d(phi*v)/dx * dxdqsi
185 // DphivDx.Multiply(dxdqsi, GradX);
186 //
187 // /*
188 // // Amplifying
189 // TPZManVector<REAL,3> minx(3,0.),maxx(3,0.);
190 //
191 // int spacedim = coord.Rows();
192 //
193 // for (int j=0; j<spacedim; j++) {
194 // minx[j] = coord.GetVal(j,0);
195 // maxx[j] = coord.GetVal(j,0);
196 // }
197 // TPZFNMatrix<6,REAL> VecMatrix(3,2,0.);
198 // for(int i = 0; i < 4; i++) {
199 // for(int j = 0; j < spacedim; j++) {
200 // minx[j] = minx[j] < coord.GetVal(j,i) ? minx[j]:coord.GetVal(j,i);
201 // maxx[j] = maxx[j] > coord.GetVal(j,i) ? maxx[j]:coord.GetVal(j,i);
202 // }
203 // }
204 // REAL delx = 0.;
205 // for (int j=0; j<spacedim; j++) {
206 // delx = delx > (maxx[j]-minx[j]) ? delx : (maxx[j]-minx[j]);
207 // }
208 //
209 // GradX *= 1./delx;
210 // */
211 //
212 // //GradX.Print("GradX");
213 //
214 // GradX.GramSchmidt(axest,jacobian);
215 // axest.Transpose(&axes);
216 // detjac = jacobian(0,0)*jacobian(1,1) - jacobian(1,0)*jacobian(0,1);
217 //
218 // if(IsZero(detjac))
219 // {
220 // detjac = ZeroTolerance();
221 // }
222 //
223 // jacinv(0,0) = jacobian(1,1)/detjac;
224 // jacinv(1,1) = jacobian(0,0)/detjac;
225 // jacinv(0,1) = -jacobian(0,1)/detjac;
226 // jacinv(1,0) = -jacobian(1,0)/detjac;
227 //
228 // /*
229 // // Left without the amplification for learning purposes
230 // jacobian *= delx;
231 // jacinv *= 1./delx;
232 // detjac *= (delx*delx);
233 // */
234 //
235 // }
236 
237  template<class T>
238  static TPZFMatrix<T> TensorProd(TPZFMatrix<T> &vec1, TPZFMatrix<T> &vec2)
239  {
240  TPZFNMatrix<9,T> res(3,3,0.);
241  for (int i = 0; i < 3; i++) {
242  for (int j = 0; j < 3; j++) {
243  res(i,j) = vec1(i,0) * vec2(j,0);
244  }
245  }
246  return res;
247  }
248 
249  // static TPZGeoEl *CreateBCGeoEl(TPZGeoEl *gel, int side,int bc);
250 
251  // /** @brief Creates a geometric element according to the type of the father element */
252  // static TPZGeoEl *CreateGeoElement(TPZGeoMesh &mesh, MElementType type,
253  // TPZVec<int64_t>& nodeindexes,
254  // int matid,
255  // int64_t& index);
256 
257  void Read(TPZStream& buf, void* context) override {
258  GeomQuad::Read(buf,context);
259  buf.Read(&fR);
260  buf.Read(fxc);
261  }
262 
263  void Write(TPZStream &buf, int withclassid) const override{
264  GeomQuad::Write(buf, withclassid);
265  buf.Write(&fR);
266  buf.Write(fxc);
267  }
268 
269  static void InsertExampleElement(TPZGeoMesh &gmesh, int matid, TPZVec<REAL> &lowercorner, TPZVec<REAL> &size);
270 
271 
272  };
273 
274 
275 }
276 
277 #endif /* defined(__PZ__TPZQuadSphere__) */
pzgeom::TPZQuadSphere::InsertExampleElement
static void InsertExampleElement(TPZGeoMesh &gmesh, int matid, TPZVec< REAL > &lowercorner, TPZVec< REAL > &size)
Definition: TPZQuadSphere.cpp:52
pzgeom::TPZQuadSphere::fxc
TPZVec< REAL > fxc
Definition: TPZQuadSphere.h:22
pzgeom::TPZQuadSphere::TPZQuadSphere
TPZQuadSphere(const TPZQuadSphere &cp, TPZGeoMesh &)
Copy constructor.
Definition: TPZQuadSphere.h:48
pzgeom::TPZQuadSphere::TPZQuadSphere
TPZQuadSphere()
Empty constructor.
Definition: TPZQuadSphere.h:32
pzgeom::TPZQuadSphere::operator=
TPZQuadSphere & operator=(const TPZQuadSphere &cp)
Definition: TPZQuadSphere.h:52
pzgeom::TPZQuadSphere::TPZQuadSphere
TPZQuadSphere(const TPZQuadSphere &cp)
Copy constructor.
Definition: TPZQuadSphere.h:43
pzgeom::TPZQuadSphere::TPZQuadSphere
TPZQuadSphere(const TPZQuadSphere &cp, std::map< int64_t, int64_t > &gl2lcNdMap)
Constructor with node map.
Definition: TPZQuadSphere.h:37
pzgeom::TPZQuadSphere::IsLinearMapping
static bool IsLinearMapping(int side)
Definition: TPZQuadSphere.h:63
pzgeom::TPZQuadSphere::X
void X(TPZFMatrix< REAL > &coord, TPZVec< T > &loc, TPZVec< T > &result) const
Definition: TPZQuadSphere.h:90
pzgeom::TPZQuadSphere::Write
void Write(TPZStream &buf, int withclassid) const override
Definition: TPZQuadSphere.h:263
pzgeom::TPZQuadSphere::Read
void Read(TPZStream &buf, void *context) override
Creates a geometric element according to the type of the father element.
Definition: TPZQuadSphere.h:257
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
pzgeom::TPZQuadSphere::TensorProd
static TPZFMatrix< T > TensorProd(TPZFMatrix< T > &vec1, TPZFMatrix< T > &vec2)
Definition: TPZQuadSphere.h:238
pzgeom::TPZQuadSphere::GradX
void GradX(TPZFMatrix< REAL > &coord, TPZVec< T > &param, TPZFMatrix< T > &gradx) const
Definition: TPZQuadSphere.h:111
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
test.res
string res
Definition: test.py:151
pzgeom::TPZQuadSphere::TypeName
static std::string TypeName()
Returns the type name of the element.
Definition: TPZQuadSphere.h:70
TPZGeoMesh
This class implements a geometric mesh for the pz environment. Geometry.
Definition: pzgmesh.h:48
TPZMatrix::Multiply
virtual void Multiply(const TPZFMatrix< TVar > &A, TPZFMatrix< TVar > &res, int opt=0) const
It mutiplies itself by TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:1916
pzgeom::TPZQuadSphere::IsGeoBlendEl
bool IsGeoBlendEl() const
declare geometry as blended element
Definition: TPZQuadSphere.cpp:101
pztopology::TPZQuadrilateral
Defines the topology of a quadrilateral element. Topology Sides 0 to 3 are vertices, sides 4 to 7 are lines, side 8 is the quadrilateral.
Definition: tpzquadrilateral.h:37
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
pzgeom::TPZQuadSphere::Top
pztopology::TPZQuadrilateral Top
Definition: TPZQuadSphere.h:25
pzgeom
Groups all classes which model the geometry.
Definition: pzgeopoint.cpp:18
pzgeom::TPZQuadSphere::TPZQuadSphere
TPZQuadSphere(TPZVec< int64_t > &nodeindexes)
Constructor with list of nodes.
Definition: TPZQuadSphere.h:27
pzgeom::TPZQuadSphere::SetData
void SetData(const REAL &R, const TPZVec< REAL > &xc)
Definition: TPZQuadSphere.h:72
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
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
Generated on Sun Aug 16 2020 11:17:10 for NeoPZ by   doxygen 1.8.13