circular_shell_mesh.template.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_CIRCULAR_SHELL_MESH_TEMPLATE_CC
31 #define OOMPH_CIRCULAR_SHELL_MESH_TEMPLATE_CC
32 
33 #include<float.h>
34 
37 
38 
39 
40 namespace oomph
41 {
42 
43 //=======================================================================
44 /// Mesh build fct
45 //=======================================================================
46  template <class ELEMENT>
48  const unsigned &nx,
49  const unsigned &ny,
50  const double &lx,
51  const double &ly)
52  {
53 
54  // Mesh can only be built with 2D Qelements.
55  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
56 
57  //Find out how many nodes there are
58  unsigned n_node = nnode();
59 
60  //Now in this case it is the Lagrangian coordinates that we want to set,
61  //so we have to loop over all nodes and set them to the Eulerian
62  //coordinates that are set by the generic mesh generator
63  for(unsigned i=0;i<n_node;i++)
64  {
65  node_pt(i)->xi(0) = scaled_x(node_pt(i)->x(0));
66  node_pt(i)->xi(1) = node_pt(i)->x(1);
67  }
68 
69 
70  // Loop over elements and nodes to find out min axial spacing
71  // for all nodes
72 
73  // Initialise map
74  std::map<Node*,double> min_dx;
75  unsigned nnod=nnode();
76  for (unsigned j=0;j<nnod;j++) min_dx[node_pt(j)]=DBL_MAX;
77 
78  // Loop over elements
79  unsigned nelem=nelement();
80  for (unsigned e=0;e<nelem;e++)
81  {
82  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(element_pt(e));
83  unsigned n_node=el_pt->nnode();
84  for(unsigned j=0;j<n_node;j++)
85  {
86  SolidNode* nod_pt=dynamic_cast<SolidNode*>(el_pt->node_pt(j));
87  double x=nod_pt->xi(0);
88  for(unsigned k=0;k<n_node;k++)
89  {
90  double dx=fabs(x-dynamic_cast<SolidNode*>(el_pt->node_pt(k))->xi(0));
91  if (dx<min_dx[nod_pt])
92  {
93  if (dx!=0.0) min_dx[nod_pt]=dx;
94  }
95  }
96  }
97  }
98 
99 
100  //Assign gradients, etc for the Lagrangian coordinates of
101  //Hermite-type elements
102 
103  //Read out number of position dofs
104  unsigned n_position_type = finite_element_pt(0)->nnodal_position_type();
105 
106  // Assign generalised Lagrangian positions (min slope, c.f. M. Heil's PhD
107  // thesis
108  if(n_position_type > 1)
109  {
110  // Default spacing
111  double xstep = (this->Xmax - this->Xmin)/((this->Np-1)*this->Nx);
112  double ystep = (this->Ymax - this->Ymin)/((this->Np-1)*this->Ny);
113 
114  // Adjust for non-uniform spacing
115  for(unsigned j=0;j<n_node;j++)
116  {
117  SolidNode* nod_pt=node_pt(j);
118 
119  // Get min. spacing for non-uniform axial spacing
120  xstep=min_dx[nod_pt];
121 
122  //The factor 0.5 is because our reference element has length 2.0
123  nod_pt->xi_gen(1,0) = 0.5*xstep;
124  nod_pt->xi_gen(2,1) = 0.5*ystep;
125  }
126  }
127 
128 
129  }
130 
131 
132 //=======================================================================
133 /// Set the undeformed coordinates of the nodes
134 //=======================================================================
135 template <class ELEMENT>
137 assign_undeformed_positions(GeomObject* const &undeformed_midplane_pt)
138 {
139 
140  // Loop over nodes in elements
141  unsigned nelem=nelement();
142  for (unsigned e=0;e<nelem;e++)
143  {
144  ELEMENT* el_pt=dynamic_cast<ELEMENT*>(element_pt(e));
145  unsigned n_node=el_pt->nnode();
146  for(unsigned n=0;n<n_node;n++)
147  {
148  //Get the Lagrangian coordinates
149  Vector<double> xi(2);
150  xi[0] = dynamic_cast<SolidNode*>(el_pt->node_pt(n))->xi(0);
151  xi[1] = dynamic_cast<SolidNode*>(el_pt->node_pt(n))->xi(1);
152 
153  //Assign memory for values of derivatives, etc
154  Vector<double> R(3);
155  DenseMatrix<double> a(2,3);
156  RankThreeTensor<double> dadxi(2,2,3);
157 
158  //Get the geometrical information from the geometric object
159  undeformed_midplane_pt->d2position(xi,R,a,dadxi);
160 
161 
162  // Get derivatives of Lagr coordinates w.r.t. local coords
163  DenseMatrix<double> dxids(2,2);
164  Vector<double> s(2);
165  el_pt->local_coordinate_of_node(n,s);
166  el_pt->interpolated_dxids(s,dxids);
167  double dxds=dxids(0,0);
168 
169  //Loop over coordinate directions
170  for(unsigned i=0;i<3;i++)
171  {
172  //Set the position
173  el_pt->node_pt(n)->x_gen(0,i) = R[i];
174 
175  //Set generalised positions
176  el_pt->node_pt(n)->x_gen(1,i)=a(0,i)*dxds;
177  el_pt->node_pt(n)->x_gen(2,i)=
178  0.5*a(1,i)*((this->Ymax-this->Ymin)/this->Ny);
179 
180  // Set the mixed derivative
181  el_pt->node_pt(n)->x_gen(3,i) = 0.0;
182 
183  // Check for warping
184  if (dadxi(0,1,i)!=0.0)
185  {
186 
187  std::ostringstream error_stream;
188  error_stream
189  << "Undef. GeomObject for this shell mesh should not be warped!\n"
190  << "It may be possible to generalise the mesh generator to \n"
191  << "deal with this case -- feel free to have a go...\n";
192  throw OomphLibError(
193  error_stream.str(),
194  OOMPH_CURRENT_FUNCTION,
195  OOMPH_EXCEPTION_LOCATION);
196  }
197  }
198  }
199  }
200 }
201 
202 
203 }
204 
205 
206 
207 
208 
209 // namespace oomph
210 // {
211 
212 
213 
214 // //=======================================================================
215 // /// Mesh constructor
216 // /// Argument list:
217 // /// nx : number of elements in the axial direction
218 // /// ny : number of elements in the azimuthal direction
219 // /// lx : length in the axial direction
220 // /// ly : length in theta direction
221 // //=======================================================================
222 // template <class ELEMENT>
223 // CircularCylindricalShellMesh<ELEMENT>::CircularCylindricalShellMesh(
224 // const unsigned &nx,
225 // const unsigned &ny,
226 // const double &lx,
227 // const double &ly,
228 // TimeStepper* time_stepper_pt) :
229 // RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt)
230 // {
231 // //Find out how many nodes there are
232 // unsigned n_node = nnode();
233 
234 // //Now in this case it is the Lagrangian coordinates that we want to set,
235 // //so we have to loop over all nodes and set them to the Eulerian
236 // //coordinates that are set by the generic mesh generator
237 // for(unsigned i=0;i<n_node;i++)
238 // {
239 // node_pt(i)->xi(0) = node_pt(i)->x(0);
240 // node_pt(i)->xi(1) = node_pt(i)->x(1);
241 // }
242 
243 
244 // //Assign gradients, etc for the Lagrangian coordinates of
245 // //Hermite-type elements
246 
247 // //Read out number of position dofs
248 // unsigned n_position_type = finite_element_pt(0)->nnodal_position_type();
249 
250 // //If this is greater than 1 set the slopes, which are the distances between
251 // //nodes. If the spacing were non-uniform, this part would be more difficult
252 // if(n_position_type > 1)
253 // {
254 // double xstep = (this->Xmax - this->Xmin)/((this->Np-1)*this->Nx);
255 // double ystep = (this->Ymax - this->Ymin)/((this->Np-1)*this->Ny);
256 // for(unsigned n=0;n<n_node;n++)
257 // {
258 // //The factor 0.5 is because our reference element has length 2.0
259 // node_pt(n)->xi_gen(1,0) = 0.5*xstep;
260 // node_pt(n)->xi_gen(2,1) = 0.5*ystep;
261 // }
262 // }
263 // }
264 
265 
266 // //=======================================================================
267 // /// Set the undeformed coordinates of the nodes
268 // //=======================================================================
269 // template <class ELEMENT>
270 // void CircularCylindricalShellMesh<ELEMENT>::assign_undeformed_positions(
271 // GeomObject* const &undeformed_midplane_pt)
272 // {
273 // //Find out how many nodes there are
274 // unsigned n_node = nnode();
275 
276 // //Loop over all the nodes
277 // for(unsigned n=0;n<n_node;n++)
278 // {
279 // //Get the Lagrangian coordinates
280 // Vector<double> xi(2);
281 // xi[0] = node_pt(n)->xi(0);
282 // xi[1] = node_pt(n)->xi(1);
283 
284 // //Assign memory for values of derivatives, etc
285 // Vector<double> R(3);
286 // DenseMatrix<double> a(2,3);
287 // RankThreeTensor<double> dadxi(2,2,3);
288 
289 // //Get the geometrical information from the geometric object
290 // undeformed_midplane_pt->d2position(xi,R,a,dadxi);
291 
292 // //Loop over coordinate directions
293 // for(unsigned i=0;i<3;i++)
294 // {
295 // //Set the position
296 // node_pt(n)->x_gen(0,i) = R[i];
297 
298 // //Set the derivative wrt Lagrangian coordinates
299 // //Note that we need to scale by the length of each element here!!
300 // node_pt(n)->x_gen(1,i) = 0.5*a(0,i)*((this->Xmax - this->Xmin)/this->Nx);
301 // node_pt(n)->x_gen(2,i) = 0.5*a(1,i)*((this->Ymax - this->Ymin)/this->Ny);
302 
303 // //Set the mixed derivative
304 // //(symmetric so doesn't matter which one we use)
305 // node_pt(n)->x_gen(3,i) = 0.25*dadxi(0,1,i);
306 // }
307 // }
308 // }
309 
310 
311 // }
312 #endif
cstr elem_len * i
Definition: cfortran.h:607
e
Definition: cfortran.h:575
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
Definition: nodes.h:1740
double & xi_gen(const unsigned &k, const unsigned &i)
Reference to the generalised Lagrangian position. `Type': k; 'Coordinate direction: i...
Definition: nodes.h:1759
A Rank 3 Tensor class.
Definition: matrices.h:1337
static char t char * s
Definition: cfortran.h:572
void build_mesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly)
Mesh build helper fct.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
virtual void d2position(const Vector< double > &zeta, RankThreeTensor< double > &ddrdzeta) const
2nd derivative of position Vector w.r.t. to coordinates: = ddrdzeta(alpha,beta,i). Evaluated at current time.
Definition: geom_objects.h:313
void assign_undeformed_positions(GeomObject *const &undeformed_midplane_pt)
In all elastic problems, the nodes must be assigned an undeformed, or reference, position, corresponding to the stress-free state of the elastic body. This function assigns the undeformed position for the nodes on the elastic tube.