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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26
27// Include guard to prevent multiple inclusions of the header
28#ifndef OOMPH_FACE_MESH_PROJECT_HEADER
29#define OOMPH_FACE_MESH_PROJECT_HEADER
30
31namespace oomph
32{
33 /// ///////////////////////////////////////////////////////////////////////
34 /// ///////////////////////////////////////////////////////////////////////
35 /// ///////////////////////////////////////////////////////////////////////
36
37
38 //======================================================================
39 /// Class that makes the finite element specified as template argument
40 /// projectable -- on the assumption that all fields are interpolated
41 /// by isoparametric Lagrange interpolation between the nodes.
42 //======================================================================
43 template<class ELEMENT>
45 : public virtual ProjectableElement<ELEMENT>
46 {
47 public:
48 /// Constructor
50 {
51 Boundary_id = UINT_MAX;
52 }
53
54 /// Nodal value of boundary coordinate
55 double zeta_nodal(const unsigned& n,
56 const unsigned& k,
57 const unsigned& i) const
58 {
59 // Vector in which to hold the intrinsic coordinate
60 Vector<double> zeta(this->dim());
61
62#ifdef PARANOID
63 if (Boundary_id == UINT_MAX)
64 {
65 std::ostringstream error_message;
66 error_message << "Boundary_id is (still) UINT_MAX -- please set\n"
67 << "the actual value with set_boundary_id(...)\n";
68 throw OomphLibError(error_message.str(),
69 OOMPH_CURRENT_FUNCTION,
70 OOMPH_EXCEPTION_LOCATION);
71 }
72#endif
73
74 // Get the k-th generalised boundary coordinate at node n
76
77 // Return the individual coordinate
78 return zeta[i];
79 }
80
81 /// Boundary id
82 void set_boundary_id(const unsigned& boundary_id)
83 {
85 }
86
87 /// Boundary id
88 unsigned boundary_id() const
89 {
90#ifdef PARANOID
91 if (Boundary_id == UINT_MAX)
92 {
93 std::ostringstream error_message;
94 error_message << "Boundary_id is (still) UINT_MAX -- please set\n"
95 << "the actual value with set_boundary_id(...)\n";
96 throw OomphLibError(error_message.str(),
97 OOMPH_CURRENT_FUNCTION,
98 OOMPH_EXCEPTION_LOCATION);
99 }
100#endif
101 return Boundary_id;
102 }
103
104 /// Specify the values associated with field fld.
105 /// The information is returned in a vector of pairs which comprise
106 /// the Data object and the value within it, that correspond to field fld.
108 {
109 // Create the vector
110 unsigned nnod = this->nnode();
111 Vector<std::pair<Data*, unsigned>> data_values(nnod);
112
113 // Loop over all nodes
114 for (unsigned j = 0; j < nnod; j++)
115 {
116 // Add the data value associated field: The node itself
117 data_values[j] = std::make_pair(this->node_pt(j), fld);
118 }
119
120 // Return the vector
121 return data_values;
122 }
123
124 /// Number of fields to be projected.
126 {
127 return this->node_pt(0)->nvalue();
128 }
129
130 /// Number of history values to be stored for fld-th field
131 /// (includes current value!). Extract from first node but assume it's
132 /// the same for all.
133 unsigned nhistory_values_for_projection(const unsigned& fld)
134 {
135 return this->node_pt(0)->ntstorage();
136 }
137
138 /// Number of positional history values (Note: count includes
139 /// current value!). Extract from first node but assume it's
140 /// the same for all.
142 {
143 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
144 }
145
146
147 /// Return Jacobian of mapping and shape functions of field fld
148 /// at local coordinate s.
149 double jacobian_and_shape_of_field(const unsigned& fld,
150 const Vector<double>& s,
151 Shape& psi)
152 {
153 this->shape(s, psi);
154 return this->J_eulerian(s);
155 }
156
157
158 /// Return interpolated field fld at local coordinate s, at time
159 /// level t (t=0: present; t>0: history values)
160 double get_field(const unsigned& t,
161 const unsigned& fld,
162 const Vector<double>& s)
163 {
164 // Local shape function
165 unsigned n_node = this->nnode();
166 Shape psi(n_node);
167
168 // Find values of shape function
169 this->shape(s, psi);
170
171 // Initialise value of u
172 double interpolated_u = 0.0;
173
174 // Sum over the local nodes
175 for (unsigned l = 0; l < n_node; l++)
176 {
177 interpolated_u += this->nodal_value(t, l, fld) * psi[l];
178 }
179 return interpolated_u;
180 }
181
182 /// Return number of values in field fld
183 unsigned nvalue_of_field(const unsigned& fld)
184 {
185 return this->nnode();
186 }
187
188
189 /// Return local equation number of value j in field fld. Assumed to
190 /// be the local nodal equation.
191 int local_equation(const unsigned& fld, const unsigned& j)
192 {
193 return this->nodal_local_eqn(j, fld);
194 }
195
196
197 private:
198 /// Boundary id
199 unsigned Boundary_id;
200 };
201
202
203 /// ////////////////////////////////////////////////////////////////////
204 /// ////////////////////////////////////////////////////////////////////
205 /// ////////////////////////////////////////////////////////////////////
206
207
208 //======================================================================
209 /// Class that makes backup (via a deep copy) of a mesh, keeping alive
210 /// enough information to allow the solution that is currently stored
211 /// on the mesh to be projected onto another mesh sometime in the
212 /// future (when the original mesh may already have been deleted).
213 /// This is mainly useful for the projection of additional
214 /// nodal values (such as Lagrange multipliers) created by FaceElements.
215 /// ASSUMPTION: All fields in the element are represented by isoparametric
216 /// Lagrange interpolation between the nodal values.
217 /// Any fields that do not fall into this category will not
218 /// be copied across correctly and if you're unlucky
219 /// the code may die...).
220 //======================================================================
221 template<class GEOMETRIC_ELEMENT>
222 class BackupMeshForProjection : public virtual Mesh
223 {
224 public:
225 /// Constructor: Pass existing mesh and the boundary ID (need to find
226 /// the boundary coordinates that are used for the projection.
227 /// Optional final argument specifies the ID of the field (i.e. the
228 /// index of the relevant nodal value!) to be projected.
229 /// If omitted, we project all of them.
231 Mesh* mesh_pt,
232 const unsigned& boundary_id,
233 const unsigned& id_of_field_to_be_projected = UINT_MAX)
234 : Boundary_id(boundary_id),
235 ID_of_field_to_be_projected(id_of_field_to_be_projected)
236 {
237 // Find unique nodes (via elements because Node_pt vector in
238 // original mesh may not have been filled (typical for most
239 // face element meshes)
240 unsigned nel = mesh_pt->nelement();
241 Element_pt.reserve(nel);
242 Node_pt.reserve(mesh_pt->nnode());
243 for (unsigned e = 0; e < nel; e++)
244 {
245 FiniteElement* el_pt = mesh_pt->finite_element_pt(e);
246 if (el_pt != 0)
247 {
248 // Make new element
249#ifdef PARANOID
250 if (dynamic_cast<GEOMETRIC_ELEMENT*>(mesh_pt->element_pt(e)) == 0)
251 {
252 std::ostringstream error_message;
253 error_message << "Element is of wrong type " << typeid(el_pt).name()
254 << " doesn't match template parameter!" << std::endl;
255 throw OomphLibError(error_message.str(),
256 OOMPH_CURRENT_FUNCTION,
257 OOMPH_EXCEPTION_LOCATION);
258
259 if (el_pt->ninternal_data() != 0)
260 {
261 std::ostringstream error_message;
262 error_message
263 << "Internal data will NOT be projected across!\n"
264 << "If you want this functionality you'll have to \n"
265 << "implement it yourself" << std::endl;
266 OomphLibWarning(error_message.str(),
267 OOMPH_CURRENT_FUNCTION,
268 OOMPH_EXCEPTION_LOCATION);
269 }
270 }
271#endif
272
273 // Make a new element
276 GEOMETRIC_ELEMENT>;
277
278 // Set boundary ID
279 new_el_pt->set_boundary_id(Boundary_id);
280
281 // Set nodal dimension
282 unsigned nodal_dim = el_pt->node_pt(0)->ndim();
283 new_el_pt->set_nodal_dimension(nodal_dim);
284
285 // Add it to mesh
286 add_element_pt(new_el_pt);
287
288 // Create new nodes if needed
289 unsigned nnod = el_pt->nnode();
290 for (unsigned j = 0; j < nnod; j++)
291 {
292 Node* old_node_pt = el_pt->node_pt(j);
293 if (New_node_pt[old_node_pt] == 0)
294 {
295 Node* new_nod_pt = 0;
296
297
298#ifdef PARANOID
299 // Check boundary node-ness
300 if (!old_node_pt->is_on_boundary())
301 {
302 std::ostringstream error_message;
303 error_message << "Node isn't on a boundary!" << std::endl;
304 throw OomphLibError(error_message.str(),
305 OOMPH_CURRENT_FUNCTION,
306 OOMPH_EXCEPTION_LOCATION);
307 }
308#endif
309
310
311 // How many values are we projecting? Default: All
312 unsigned nval = old_node_pt->nvalue();
313 unsigned first_index_in_old_node = 0;
314 if (ID_of_field_to_be_projected != UINT_MAX)
315 {
316 nval = dynamic_cast<BoundaryNodeBase*>(old_node_pt)
317 ->nvalue_assigned_by_face_element(
319 first_index_in_old_node =
320 dynamic_cast<BoundaryNodeBase*>(old_node_pt)
321 ->index_of_first_value_assigned_by_face_element(
323 }
324
325 // Build new node
326 new_nod_pt =
327 new BoundaryNode<Node>(old_node_pt->time_stepper_pt(),
328 old_node_pt->ndim(),
329 old_node_pt->nposition_type(),
330 nval);
331
332 // Copy data across
333 unsigned n_time = old_node_pt->ntstorage();
334 for (unsigned t = 0; t < n_time; t++)
335 {
336 for (unsigned i = 0; i < nval; i++)
337 {
338 new_nod_pt->set_value(
339 t, i, old_node_pt->value(t, first_index_in_old_node + i));
340 }
341 }
342
343 // Copy nodal positions
344 unsigned n_dim = old_node_pt->ndim();
345 for (unsigned i = 0; i < n_dim; i++)
346 {
347 new_nod_pt->x(i) = old_node_pt->x(i);
348 }
349
350 // Add to boundary
351 new_nod_pt->add_to_boundary(Boundary_id);
352
353 // Transfer boundary coordinates
354#ifdef PARANOID
355 if (!old_node_pt->is_on_boundary(Boundary_id))
356 {
357 std::ostringstream error_message;
358 error_message << "Boundary ID specified as " << Boundary_id
359 << " but node isn't actually on that boundary!"
360 << std::endl;
361 throw OomphLibError(error_message.str(),
362 OOMPH_CURRENT_FUNCTION,
363 OOMPH_EXCEPTION_LOCATION);
364 }
365#endif
366
367 // Get vector of coordinates on mesh boundary from old node
368 unsigned n = old_node_pt->ncoordinates_on_boundary(Boundary_id);
369 Vector<double> zeta(n);
370 old_node_pt->get_coordinates_on_boundary(Boundary_id, zeta);
371
372 // Set for new node
373 new_nod_pt->set_coordinates_on_boundary(Boundary_id, zeta);
374
375 // Add node
376 Node_pt.push_back(new_nod_pt);
377
378 // Setup association
379 New_node_pt[old_node_pt] = new_nod_pt;
380 }
381
382 // Set node pointer from new element
383 new_el_pt->node_pt(j) = New_node_pt[old_node_pt];
384 }
385 }
386 }
387 }
388
389 /// Project the solution that was present in the original mesh
390 /// and from which this backup mesh was created onto the mesh
391 /// pointed to by new_mesh_pt. Note that elements in the new mesh do
392 /// not have to be projectable. The original mesh may by now have
393 /// been deleted.
394 void project_onto_new_mesh(Mesh* new_mesh_pt)
395 {
396 // Make copy of new mesh that we can project onto
397 BackupMeshForProjection<GEOMETRIC_ELEMENT>* projectable_new_mesh_pt =
400
401 // Create projection problem
404 proj_problem_pt = new ProjectionProblem<
406
407 // Set the mesh we want to project onto
408 proj_problem_pt->mesh_pt() = projectable_new_mesh_pt;
409
410 // Project from projectable copy of original mesh -- this one!
411 bool dont_project_positions = true;
412 proj_problem_pt->project(this, dont_project_positions);
413
414 // Copy nodal values onto the corresponding nodes in the new mesh
415 projectable_new_mesh_pt->copy_onto_original_mesh();
416
417 // Kill!
418 delete proj_problem_pt;
419 delete projectable_new_mesh_pt;
420 }
421
422
423 /// Copy nodal values back onto original mesh from which this
424 /// mesh was built. This obviously only makes sense if the original
425 /// mesh still exists!
427 {
428 for (std::map<Node*, Node*>::iterator it = New_node_pt.begin();
429 it != New_node_pt.end();
430 it++)
431 {
432 // Get old node (in the previously existing mesh)
433 Node* old_node_pt = (*it).first;
434
435 // ...and corresponding new one (in the mesh where we did the projection
436 Node* new_node_pt = (*it).second;
437
438 // How many values are we moving across?
439 unsigned nval = old_node_pt->nvalue();
440 unsigned first_index_in_old_node = 0;
441 if (ID_of_field_to_be_projected != UINT_MAX)
442 {
443 nval =
444 dynamic_cast<BoundaryNodeBase*>(old_node_pt)
445 ->nvalue_assigned_by_face_element(ID_of_field_to_be_projected);
446 first_index_in_old_node =
447 dynamic_cast<BoundaryNodeBase*>(old_node_pt)
448 ->index_of_first_value_assigned_by_face_element(
450 }
451
452 // Copy data across (include offset in orig mesh)
453 unsigned n_time = old_node_pt->ntstorage();
454 for (unsigned t = 0; t < n_time; t++)
455 {
456 for (unsigned i = 0; i < nval; i++)
457 {
458 old_node_pt->set_value(
459 t, first_index_in_old_node + i, new_node_pt->value(t, i));
460 }
461 }
462 }
463 }
464
465
466 private:
467 /// Map returning new node, labeled by node point in original mesh
468 std::map<Node*, Node*> New_node_pt;
469
470 /// Boundary id
471 unsigned Boundary_id;
472
473 /// ID of field to be projected (UINT_MAX if there isn't one, in
474 /// which case we're doing all of them.
476 };
477
478
479 /// //////////////////////////////////////////////////////////////////
480 /// //////////////////////////////////////////////////////////////////
481 /// //////////////////////////////////////////////////////////////////
482
483} // namespace oomph
484
485#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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...
std::map< Node *, Node * > New_node_pt
Map returning new node, labeled by node point in original mesh.
void copy_onto_original_mesh()
Copy nodal values back onto original mesh from which this mesh was built. This obviously only makes s...
unsigned Boundary_id
Boundary id.
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.
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...
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition: nodes.h:1996
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh....
Definition: nodes.h:2242
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:879
A general Finite Element class.
Definition: elements.h:1313
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:4103
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:1390
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:2593
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...
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:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:823
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!). Extract from first node bu...
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...
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
Nodal value of boundary coordinate.
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 ...
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.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
unsigned nfields_for_projection()
Number of fields to be projected.
void set_boundary_id(const unsigned &boundary_id)
Boundary id.
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.
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 general mesh class.
Definition: mesh.h:67
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:448
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:617
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
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:2363
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
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:2336
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:2394
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
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....
Definition: nodes.cc:2379
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
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:2408
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
Projection problem. This is created during the adaptation of unstructured meshes and it is assumed th...
Definition: projection.h:695
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
//////////////////////////////////////////////////////////////////// ////////////////////////////////...