face_mesh_project.h
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: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
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 
31 //Include guard to prevent multiple inclusions of the header
32 #ifndef OOMPH_FACE_MESH_PROJECT_HEADER
33 #define OOMPH_FACE_MESH_PROJECT_HEADER
34 
35 namespace oomph
36 {
37 
38 
39 //////////////////////////////////////////////////////////////////////////
40 //////////////////////////////////////////////////////////////////////////
41 //////////////////////////////////////////////////////////////////////////
42 
43 
44  //======================================================================
45  /// \short Class that makes the finite element specified as template argument
46  /// projectable -- on the assumption that all fields are interpolated
47  /// by isoparametric Lagrange interpolation between the nodes.
48  //======================================================================
49  template <class ELEMENT>
51  public virtual ProjectableElement<ELEMENT>
52  {
53 
54  public:
55 
56  /// Constructor
58  {
59  Boundary_id=UINT_MAX;
60  }
61 
62  /// \short Nodal value of boundary coordinate
63  double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i)
64  const
65  {
66  //Vector in which to hold the intrinsic coordinate
67  Vector<double> zeta(this->dim());
68 
69 #ifdef PARANOID
70  if (Boundary_id==UINT_MAX)
71  {
72  std::ostringstream error_message;
73  error_message
74  << "Boundary_id is (still) UINT_MAX -- please set\n"
75  << "the actual value with set_boundary_id(...)\n";
76  throw OomphLibError(error_message.str(),
77  OOMPH_CURRENT_FUNCTION,
78  OOMPH_EXCEPTION_LOCATION);
79  }
80 #endif
81 
82  //Get the k-th generalised boundary coordinate at node n
84  Boundary_id,k,zeta);
85 
86  //Return the individual coordinate
87  return zeta[i];
88  }
89 
90  /// Boundary id
91  void set_boundary_id(const unsigned& boundary_id)
92  {
94  }
95 
96  /// Boundary id
97  unsigned boundary_id() const
98  {
99 #ifdef PARANOID
100  if (Boundary_id==UINT_MAX)
101  {
102  std::ostringstream error_message;
103  error_message
104  << "Boundary_id is (still) UINT_MAX -- please set\n"
105  << "the actual value with set_boundary_id(...)\n";
106  throw OomphLibError(error_message.str(),
107  OOMPH_CURRENT_FUNCTION,
108  OOMPH_EXCEPTION_LOCATION);
109  }
110 #endif
111  return Boundary_id;
112  }
113 
114  /// \short Specify the values associated with field fld.
115  /// The information is returned in a vector of pairs which comprise
116  /// the Data object and the value within it, that correspond to field fld.
118  {
119  // Create the vector
120  unsigned nnod=this->nnode();
121  Vector<std::pair<Data*,unsigned> > data_values(nnod);
122 
123  // Loop over all nodes
124  for (unsigned j=0;j<nnod;j++)
125  {
126  // Add the data value associated field: The node itself
127  data_values[j]=std::make_pair(this->node_pt(j),fld);
128  }
129 
130  // Return the vector
131  return data_values;
132  }
133 
134  /// \short Number of fields to be projected.
136  {
137  return this->node_pt(0)->nvalue();
138  }
139 
140  /// \short Number of history values to be stored for fld-th field
141  /// (includes current value!). Extract from first node but assume it's
142  /// the same for all.
143  unsigned nhistory_values_for_projection(const unsigned &fld)
144  {
145  return this->node_pt(0)->ntstorage();
146  }
147 
148  ///\short Number of positional history values (Note: count includes
149  /// current value!). Extract from first node but assume it's
150  /// the same for all.
152  {
153  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
154  }
155 
156 
157  /// \short Return Jacobian of mapping and shape functions of field fld
158  /// at local coordinate s.
159  double jacobian_and_shape_of_field(const unsigned &fld,
160  const Vector<double> &s,
161  Shape &psi)
162  {
163  this->shape(s,psi);
164  return this->J_eulerian(s);
165  }
166 
167 
168 
169  /// \short Return interpolated field fld at local coordinate s, at time level
170  /// t (t=0: present; t>0: history values)
171  double get_field(const unsigned &t,
172  const unsigned &fld,
173  const Vector<double>& s)
174  {
175  //Local shape function
176  unsigned n_node=this->nnode();
177  Shape psi(n_node);
178 
179  //Find values of shape function
180  this->shape(s,psi);
181 
182  //Initialise value of u
183  double interpolated_u = 0.0;
184 
185  //Sum over the local nodes
186  for(unsigned l=0;l<n_node;l++)
187  {
188  interpolated_u += this->nodal_value(t,l,fld)*psi[l];
189  }
190  return interpolated_u;
191  }
192 
193  /// \short Return number of values in field fld
194  unsigned nvalue_of_field(const unsigned &fld)
195  {
196  return this->nnode();
197  }
198 
199 
200  /// \short Return local equation number of value j in field fld. Assumed to be
201  /// the local nodal equation.
202  int local_equation(const unsigned &fld,
203  const unsigned &j)
204  {
205  return this->nodal_local_eqn(j,fld);
206  }
207 
208 
209  private:
210 
211  /// Boundary id
212  unsigned Boundary_id;
213 
214  };
215 
216 
217 
218 ///////////////////////////////////////////////////////////////////////
219 ///////////////////////////////////////////////////////////////////////
220 ///////////////////////////////////////////////////////////////////////
221 
222 
223  //======================================================================
224  /// Class that makes backup (via a deep copy) of a mesh, keeping alive
225  /// enough information to allow the solution that is currently stored
226  /// on the mesh to be projected onto another mesh sometime in the
227  /// future (when the original mesh may already have been deleted).
228  /// This is mainly useful for the projection of additional
229  /// nodal values (such as Lagrange multipliers) created by FaceElements.
230  /// ASSUMPTION: All fields in the element are represented by isoparametric
231  /// Lagrange interpolation between the nodal values.
232  /// Any fields that do not fall into this category will not
233  /// be copied across correctly and if you're unlucky
234  /// the code may die...).
235  //======================================================================
236  template <class GEOMETRIC_ELEMENT>
237  class BackupMeshForProjection : public virtual Mesh
238  {
239 
240  public:
241 
242  /// \short Constructor: Pass existing mesh and the boundary ID (need to find
243  /// the boundary coordinates that are used for the projection.
244  /// Optional final argument specifies the ID of the field (i.e. the
245  /// index of the relevant nodal value!) to be projected.
246  /// If omitted, we project all of them.
247  BackupMeshForProjection(Mesh* mesh_pt, const unsigned& boundary_id,
248  const unsigned& id_of_field_to_be_projected=UINT_MAX)
249  : Boundary_id(boundary_id),
250  ID_of_field_to_be_projected(id_of_field_to_be_projected)
251  {
252  // Find unique nodes (via elements because Node_pt vector in
253  // original mesh may not have been filled (typical for most
254  // face element meshes)
255  unsigned nel=mesh_pt->nelement();
256  Element_pt.reserve(nel);
257  Node_pt.reserve(mesh_pt->nnode());
258  for (unsigned e=0;e<nel;e++)
259  {
260  FiniteElement* el_pt=mesh_pt->finite_element_pt(e);
261  if (el_pt!=0)
262  {
263  // Make new element
264 #ifdef PARANOID
265  if (dynamic_cast<GEOMETRIC_ELEMENT*>(mesh_pt->element_pt(e))==0)
266  {
267  std::ostringstream error_message;
268  error_message
269  << "Element is of wrong type " << typeid(el_pt).name()
270  << " doesn't match template parameter!" << std::endl;
271  throw OomphLibError(error_message.str(),
272  OOMPH_CURRENT_FUNCTION,
273  OOMPH_EXCEPTION_LOCATION);
274 
275  if (el_pt->ninternal_data()!=0)
276  {
277  std::ostringstream error_message;
278  error_message
279  << "Internal data will NOT be projected across!\n"
280  << "If you want this functionality you'll have to \n"
281  << "implement it yourself" << std::endl;
282  OomphLibWarning(error_message.str(),
283  OOMPH_CURRENT_FUNCTION,
284  OOMPH_EXCEPTION_LOCATION);
285  }
286  }
287 #endif
288 
289  // Make a new element
292 
293  // Set boundary ID
294  new_el_pt->set_boundary_id(Boundary_id);
295 
296  // Set nodal dimension
297  unsigned nodal_dim=el_pt->node_pt(0)->ndim();
298  new_el_pt->set_nodal_dimension(nodal_dim);
299 
300  // Add it to mesh
301  add_element_pt(new_el_pt);
302 
303  // Create new nodes if needed
304  unsigned nnod=el_pt->nnode();
305  for (unsigned j=0;j<nnod;j++)
306  {
307  Node* old_node_pt=el_pt->node_pt(j);
308  if (New_node_pt[old_node_pt]==0)
309  {
310  Node* new_nod_pt=0;
311 
312 
313 #ifdef PARANOID
314  // Check boundary node-ness
315  if (!old_node_pt->is_on_boundary())
316  {
317  std::ostringstream error_message;
318  error_message
319  << "Node isn't on a boundary!"
320  << std::endl;
321  throw OomphLibError(error_message.str(),
322  OOMPH_CURRENT_FUNCTION,
323  OOMPH_EXCEPTION_LOCATION);
324  }
325 #endif
326 
327 
328  // How many values are we projecting? Default: All
329  unsigned nval=old_node_pt->nvalue();
330  unsigned first_index_in_old_node=0;
331  if (ID_of_field_to_be_projected!=UINT_MAX)
332  {
333  nval=dynamic_cast<BoundaryNodeBase*>
334  (old_node_pt)->nvalue_assigned_by_face_element
336  first_index_in_old_node=dynamic_cast<BoundaryNodeBase*>
337  (old_node_pt)->index_of_first_value_assigned_by_face_element
339  }
340 
341  // Build new node
342  new_nod_pt=new BoundaryNode<Node>(old_node_pt->time_stepper_pt(),
343  old_node_pt->ndim(),
344  old_node_pt->nposition_type(),
345  nval);
346 
347  // Copy data across
348  unsigned n_time = old_node_pt->ntstorage();
349  for(unsigned t=0;t<n_time;t++)
350  {
351  for(unsigned i=0;i<nval;i++)
352  {
353  new_nod_pt->set_value
354  (t,i,old_node_pt->value(t,first_index_in_old_node+i));
355  }
356  }
357 
358  // Copy nodal positions
359  unsigned n_dim=old_node_pt->ndim();
360  for(unsigned i=0;i<n_dim;i++)
361  {
362  new_nod_pt->x(i)=old_node_pt->x(i);
363  }
364 
365  // Add to boundary
366  new_nod_pt->add_to_boundary(Boundary_id);
367 
368  // Transfer boundary coordinates
369 #ifdef PARANOID
370  if (!old_node_pt->is_on_boundary(Boundary_id))
371  {
372  std::ostringstream error_message;
373  error_message
374  << "Boundary ID specified as " << Boundary_id
375  << " but node isn't actually on that boundary!"
376  << std::endl;
377  throw OomphLibError(error_message.str(),
378  OOMPH_CURRENT_FUNCTION,
379  OOMPH_EXCEPTION_LOCATION);
380  }
381 #endif
382 
383  // Get vector of coordinates on mesh boundary from old node
384  unsigned n=old_node_pt->ncoordinates_on_boundary(Boundary_id);
385  Vector<double> zeta(n);
386  old_node_pt->get_coordinates_on_boundary(Boundary_id,zeta);
387 
388  // Set for new node
389  new_nod_pt->set_coordinates_on_boundary(Boundary_id,zeta);
390 
391  // Add node
392  Node_pt.push_back(new_nod_pt);
393 
394  // Setup association
395  New_node_pt[old_node_pt]=new_nod_pt;
396  }
397 
398  // Set node pointer from new element
399  new_el_pt->node_pt(j)=New_node_pt[old_node_pt];
400  }
401  }
402  }
403  }
404 
405  /// \short Project the solution that was present in the original mesh
406  /// and from which this backup mesh was created onto the mesh
407  /// pointed to by new_mesh_pt. Note that elements in the new mesh do
408  /// not have to be projectable. The original mesh may by now have
409  /// been deleted.
410  void project_onto_new_mesh(Mesh* new_mesh_pt)
411  {
412  // Make copy of new mesh that we can project onto
413  BackupMeshForProjection<GEOMETRIC_ELEMENT>* projectable_new_mesh_pt=
416 
417  // Create projection problem
419  <GEOMETRIC_ELEMENT> >*
420  proj_problem_pt=new
422  <GEOMETRIC_ELEMENT> >;
423 
424  // Set the mesh we want to project onto
425  proj_problem_pt->mesh_pt()=projectable_new_mesh_pt;
426 
427  // Project from projectable copy of original mesh -- this one!
428  bool dont_project_positions=true;
429  proj_problem_pt->project(this,dont_project_positions);
430 
431  // Copy nodal values onto the corresponding nodes in the new mesh
432  projectable_new_mesh_pt->copy_onto_original_mesh();
433 
434  // Kill!
435  delete proj_problem_pt;
436  delete projectable_new_mesh_pt;
437  }
438 
439 
440  /// \short Copy nodal values back onto original mesh from which this
441  /// mesh was built. This obviously only makes sense if the original
442  /// mesh still exists!
444  {
445 
446  for (std::map<Node*,Node*>::iterator it=New_node_pt.begin();
447  it!=New_node_pt.end();it++)
448  {
449  // Get old node (in the previously existing mesh)
450  Node* old_node_pt=(*it).first;
451 
452  // ...and corresponding new one (in the mesh where we did the projection
453  Node* new_node_pt=(*it).second;
454 
455  // How many values are we moving across?
456  unsigned nval=old_node_pt->nvalue();
457  unsigned first_index_in_old_node=0;
458  if (ID_of_field_to_be_projected!=UINT_MAX)
459  {
460  nval=dynamic_cast<BoundaryNodeBase*>(old_node_pt)->
461  nvalue_assigned_by_face_element
463  first_index_in_old_node=dynamic_cast<BoundaryNodeBase*>(old_node_pt)->
464  index_of_first_value_assigned_by_face_element
466  }
467 
468  // Copy data across (include offset in orig mesh)
469  unsigned n_time = old_node_pt->ntstorage();
470  for(unsigned t=0;t<n_time;t++)
471  {
472  for(unsigned i=0;i<nval;i++)
473  {
474  old_node_pt->set_value
475  (t,first_index_in_old_node+i,new_node_pt->value(t,i));
476  }
477  }
478 
479  }
480  }
481 
482 
483  private:
484 
485  /// Map returning new node, labeled by node point in original mesh
486  std::map<Node*,Node*> New_node_pt;
487 
488  /// Boundary id
489  unsigned Boundary_id;
490 
491  /// \short ID of field to be projected (UINT_MAX if there isn't one, in which
492  /// case we're doing all of them.
494 
495  };
496 
497 
498 /////////////////////////////////////////////////////////////////////
499 /////////////////////////////////////////////////////////////////////
500 /////////////////////////////////////////////////////////////////////
501 
502 }
503 
504 #endif
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2301
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:194
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
unsigned nfields_for_projection()
Number of fields to be projected.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
std::map< Node *, Node * > New_node_pt
Map returning new node, labeled by node point in original mesh.
unsigned Boundary_id
Boundary id.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
A general Finite Element class.
Definition: elements.h:1271
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
void copy_onto_original_mesh()
Copy nodal values back onto original mesh from which this mesh was built. This obviously only makes s...
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary, so return false. This will be overloaded by BoundaryNodes.
Definition: nodes.h:1293
char t
Definition: cfortran.h:572
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
BackupMeshForProjection(Mesh *mesh_pt, const unsigned &boundary_id, const unsigned &id_of_field_to_be_projected=UINT_MAX)
Constructor: Pass existing mesh and the boundary ID (need to find the boundary coordinates that are u...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned ID_of_field_to_be_projected
ID of field to be projected (UINT_MAX if there isn't one, in which case we're doing all of them...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
e
Definition: cfortran.h:575
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interface provides run-time error checking.
Definition: nodes.cc:2287
virtual unsigned ncoordinates_on_boundary(const unsigned &b)
Get the number of boundary coordinates on mesh boundary b. Broken virtual interface provides run-time...
Definition: nodes.cc:2271
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:966
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
A class that contains the information required by Nodes that are located on Mesh boundaries. A BoundaryNode of a particular type is obtained by combining a given Node with this class. By differentiating between Nodes and BoundaryNodes we avoid a lot of un-necessary storage in the bulk Nodes.
Definition: nodes.h:1854
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh...
Definition: nodes.h:67
static char t char * s
Definition: cfortran.h:572
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
void set_boundary_id(const unsigned &boundary_id)
Boundary id.
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:197
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!). Extract from first node bu...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field (includes current value!). Extract from first ...
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
void project_onto_new_mesh(Mesh *new_mesh_pt)
Project the solution that was present in the original mesh and from which this backup mesh was create...
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void set_nodal_dimension(const unsigned &nodal_dim)
Set the dimension of the nodes in the element. This will typically only be required when constructing...
Definition: elements.h:1348
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:605
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
virtual void add_to_boundary(const unsigned &b)
Broken interface for adding the node to the mesh boundary b Essentially here for error reporting...
Definition: nodes.cc:2243
Class that makes the finite element specified as template argument projectable – on the assumption th...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Nodal value of boundary coordinate.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
A general mesh class.
Definition: mesh.h:74
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld. Assumed to be the local nodal equation...
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...