mesh_as_geometric_object.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// Header file for a class that is used to represent a mesh
27// as a geometric object
28
29// Include guards to prevent multiple inclusion of the header
30#ifndef OOMPH_MESH_AS_GEOMETRIC_OBJECT_HEADER
31#define OOMPH_MESH_AS_GEOMETRIC_OBJECT_HEADER
32
33// Config header generated by autoconfig
34#ifdef HAVE_CONFIG_H
35#include <oomph-lib-config.h>
36#endif
37
38#include <float.h>
39#include <limits.h>
40
41// Include the geometric object header file
42#include "geom_objects.h"
43
44// Sample point container
46
47
49
50namespace oomph
51{
52 /// /////////////////////////////////////////////////////////////////////
53 /// /////////////////////////////////////////////////////////////////////
54 /// /////////////////////////////////////////////////////////////////////
55
56 //========================================================================
57 /// Helper namespace for MeshAsGeomObject -- its only function creates
58 /// SamplePointContainerParameters of the right type for the default sample
59 /// point container
60 //========================================================================
61 namespace MeshAsGeomObject_Helper
62 {
63 /// Default sample point container type
65
66 /// "Factory" for SamplePointContainerParameters of the right type as
67 /// selected by Default_sample_point_container_version
69 Mesh* mesh_pt,
70 SamplePointContainerParameters*& sample_point_container_parameters_pt);
71
72 } // namespace MeshAsGeomObject_Helper
73
74
75 /// /////////////////////////////////////////////////////////////////////
76 /// /////////////////////////////////////////////////////////////////////
77 /// /////////////////////////////////////////////////////////////////////
78
79
80 //========================================================================
81 /// This class provides a GeomObject representation of a given
82 /// finite element mesh. The Lagrangian coordinate is taken to be the
83 /// dimension of the (first) element in the mesh and the Eulerian
84 /// coordinate is taken to be the dimension of the (first) node in
85 /// the mesh. If there are no elements or nodes the appropriate dimensions
86 /// will be set to zero.
87 /// The constituent elements of the mesh must have their own
88 /// GeomObject representations, so they must be FiniteElements,
89 /// and they become sub-objects
90 /// in this compound GeomObject.
91 //========================================================================
93 {
94 private:
95 /// Helper function to actually build the thing
97 SamplePointContainerParameters* sample_point_container_parameters_pt)
98 {
99 Mesh_pt = sample_point_container_parameters_pt->mesh_pt();
100 if (dynamic_cast<RefineableBinArrayParameters*>(
101 sample_point_container_parameters_pt) != 0)
102 {
104 }
105 else if (dynamic_cast<NonRefineableBinArrayParameters*>(
106 sample_point_container_parameters_pt) != 0)
107 {
109 }
110#ifdef OOMPH_HAS_CGAL
111 else if (dynamic_cast<CGALSamplePointContainerParameters*>(
112 sample_point_container_parameters_pt) != 0)
113 {
115 }
116#endif
117 else
118 {
119 throw OomphLibError("Wrong sample_point_container_parameters_pt",
120 OOMPH_CURRENT_FUNCTION,
121 OOMPH_EXCEPTION_LOCATION);
122 }
123
124#ifdef OOMPH_HAS_MPI
125
126 // Set communicator
128
129#endif
130
131
132 // Storage for the Lagrangian and Eulerian dimension
133 int dim[2] = {0, 0};
134
135 // Set the Lagrangian dimension from the dimension of the first element
136 // if it exists (if not the Lagrangian dimension will be zero)
137 if (Mesh_pt->nelement() != 0)
138 {
139 dim[0] = Mesh_pt->finite_element_pt(0)->dim();
140 }
141
142 // Read out the Eulerian dimension from the first node, if it exists.
143 //(if not the Eulerian dimension will be zero);
144 if (Mesh_pt->nnode() != 0)
145 {
146 dim[1] = Mesh_pt->node_pt(0)->ndim();
147 }
148
149 // Need to do an Allreduce to ensure that the dimension is consistent
150 // even when no elements are assigned to a certain processor
151#ifdef OOMPH_HAS_MPI
152
153 // Only a problem if the mesh has been distributed
155 {
156 // Need a non-null communicator
157 if (Communicator_pt != 0)
158 {
159 int n_proc = Communicator_pt->nproc();
160 if (n_proc > 1)
161 {
162 int dim_reduce[2];
163 MPI_Allreduce(&dim,
164 &dim_reduce,
165 2,
166 MPI_INT,
167 MPI_MAX,
168 Communicator_pt->mpi_comm());
169
170 dim[0] = dim_reduce[0];
171 dim[1] = dim_reduce[1];
172 }
173 }
174 }
175#endif
176
177 // Set the Lagrangian and Eulerian dimensions within this geometric object
178 this->set_nlagrangian_and_ndim(static_cast<unsigned>(dim[0]),
179 static_cast<unsigned>(dim[1]));
180
181 // Create temporary storage for geometric Data (don't count
182 // Data twice!
183 std::set<Data*> tmp_geom_data;
184
185 // Copy all the elements in the mesh into local storage
186 // N.B. elements must be able to have a geometric object representation.
187 unsigned n_sub_object = Mesh_pt->nelement();
188 Sub_geom_object_pt.resize(n_sub_object);
189 for (unsigned e = 0; e < n_sub_object; e++)
190 {
191 // (Try to) cast to a finite element:
193 dynamic_cast<FiniteElement*>(Mesh_pt->element_pt(e));
194
195#ifdef PARANOID
196 if (Sub_geom_object_pt[e] == 0)
197 {
198 std::ostringstream error_message;
199 error_message << "Unable to dynamic cast element: " << std::endl
200 << "into a FiniteElement: GeomObject representation is "
201 "not possible\n";
202 throw OomphLibError(error_message.str(),
203 OOMPH_CURRENT_FUNCTION,
204 OOMPH_EXCEPTION_LOCATION);
205 }
206#endif
207
208 // Add the geometric Data of each element into set
209 unsigned ngeom_data = Sub_geom_object_pt[e]->ngeom_data();
210 for (unsigned i = 0; i < ngeom_data; i++)
211 {
212 tmp_geom_data.insert(Sub_geom_object_pt[e]->geom_data_pt(i));
213 }
214 }
215
216 // Now copy unique geom Data values across into vector
217 unsigned ngeom = tmp_geom_data.size();
218 Geom_data_pt.resize(ngeom);
219 typedef std::set<Data*>::iterator IT;
220 unsigned count = 0;
221 for (IT it = tmp_geom_data.begin(); it != tmp_geom_data.end(); it++)
222 {
223 Geom_data_pt[count] = *it;
224 count++;
225 }
226
227 // Build the right type of bin array
229 {
231
233 new RefineableBinArray(sample_point_container_parameters_pt);
234 break;
235
237
239 new NonRefineableBinArray(sample_point_container_parameters_pt);
240 break;
241
242#ifdef OOMPH_HAS_CGAL
243
245
247 new CGALSamplePointContainer(sample_point_container_parameters_pt);
248 break;
249
250#endif
251
252 default:
253
254 oomph_info << "Sample_point_container_version = "
255 << Sample_point_container_version << std::endl;
256 throw OomphLibError("Sample_point_container_version",
257 OOMPH_CURRENT_FUNCTION,
258 OOMPH_EXCEPTION_LOCATION);
259 }
260 }
261
262
263 /// Vector of pointers to Data items that affects the object's shape
265
266 /// Internal storage for the elements that constitute the object
268
269 /// Pointer to the sample point container
271
272#ifdef OOMPH_HAS_MPI
273
274 /// Communicator
276
277#endif
278
279 /// Pointer to mesh
281
282 /// Which version of the sample point container
283 /// are we using?
285
286 public:
287 /// Pointer to the sample point container
289 {
291 }
292
293 /// Return pointer to e-th finite element
295 {
296 return Sub_geom_object_pt[e];
297 }
298
299
300 /// Which sample point container is used in locate zeta? (uses enum
301 /// Sample_Point_Container_Type)
303 {
305 }
306
307 /// Number of elements in the underlying mesh
308 unsigned nelement()
309 {
310 return Sub_geom_object_pt.size();
311 }
312
313 /// Constructor
314 MeshAsGeomObject(Mesh* const& mesh_pt) : GeomObject()
315 {
316 // Create default parameters
317 SamplePointContainerParameters* sample_point_container_parameters_pt = 0;
319 mesh_pt, sample_point_container_parameters_pt);
320
321 // Build the bastard
322 build_it(sample_point_container_parameters_pt);
323 delete sample_point_container_parameters_pt;
324 }
325
326
327 /// Constructor
329 SamplePointContainerParameters* sample_point_container_parameters_pt)
330 : GeomObject()
331 {
332 build_it(sample_point_container_parameters_pt);
333 }
334
335 /// Empty Constructor
337
338 /// Destructor
340 {
342 }
343
344 /// Broken copy constructor
346
347 /// Broken assignment operator
348 void operator=(const MeshAsGeomObject&) = delete;
349
350 /// How many items of Data does the shape of the object depend on?
351 unsigned ngeom_data() const
352 {
353 return Geom_data_pt.size();
354 }
355
356 /// Return pointer to the j-th Data item that the object's
357 /// shape depends on
358 Data* geom_data_pt(const unsigned& j)
359 {
360 return Geom_data_pt[j];
361 }
362
363 /// Find the sub geometric object and local coordinate therein that
364 /// corresponds to the intrinsic coordinate zeta. If sub_geom_object_pt=0
365 /// on return from this function, none of the constituent sub-objects
366 /// contain the required coordinate. Following from the general
367 /// interface to this function in GeomObjects,
368 /// setting the optional bool argument to true means that each
369 /// time the sub-object's locate_zeta function is called, the coordinate
370 /// argument "s" is used as the initial guess. However, this doesn't
371 /// make sense here and the argument is ignored (though a warning
372 /// is issued when the code is compiled in PARANOID setting)
373 void locate_zeta(const Vector<double>& zeta,
374 GeomObject*& sub_geom_object_pt,
376 const bool& use_coordinate_as_initial_guess = false)
377 {
378#ifdef PARANOID
379 if (use_coordinate_as_initial_guess)
380 {
382 "Ignoring the use_coordinate_as_initial_guess argument.",
383 "MeshAsGeomObject::locate_zeta()",
384 OOMPH_EXCEPTION_LOCATION);
385 }
386#endif
387
388
389 // Do locate in sample point container
390 Sample_point_container_pt->locate_zeta(zeta, sub_geom_object_pt, s);
391 }
392
393 /// Return the position as a function of the intrinsic coordinate
394 /// zeta. This provides an (expensive!) default implementation in which we
395 /// loop over all the constituent sub-objects and check if they contain zeta
396 /// and then evaluate their position() function.
397 void position(const Vector<double>& zeta, Vector<double>& r) const
398 {
399 // Call position function at current timestep:
400 unsigned t = 0;
401 position(t, zeta, r);
402 }
403
404 /// Parametrised position on object: r(zeta). Evaluated at
405 /// previous timestep. t=0: current time; t>0: previous
406 /// timestep. This provides an (expensive!) default implementation in which
407 /// we loop over all the constituent sub-objects and check if they
408 /// contain zeta and then evaluate their position() function.
409 void position(const unsigned& t,
410 const Vector<double>& zeta,
411 Vector<double>& r) const
412 {
413 // Storage for the GeomObject that contains the zeta coordinate
414 // and the intrinsic coordinate within it.
415 GeomObject* sub_geom_object_pt;
416 const unsigned n_lagrangian = this->nlagrangian();
417 Vector<double> s(n_lagrangian);
418
419 // Find the sub object containing zeta, and the local intrinsic coordinate
420 // within it
421 const_cast<MeshAsGeomObject*>(this)->locate_zeta(
422 zeta, sub_geom_object_pt, s);
423 if (sub_geom_object_pt == 0)
424 {
425 std::ostringstream error_message;
426 error_message << "Cannot locate zeta ";
427 for (unsigned i = 0; i < n_lagrangian; i++)
428 {
429 error_message << zeta[i] << " ";
430 }
431 error_message << std::endl;
432 Mesh_pt->output("most_recent_mesh.dat");
433 throw OomphLibError(error_message.str(),
434 OOMPH_CURRENT_FUNCTION,
435 OOMPH_EXCEPTION_LOCATION);
436 }
437 // Call that sub-object's position function
438 sub_geom_object_pt->position(t, s, r);
439
440 } // end of position
441
442 /// Return the derivative of the position
443 void dposition(const Vector<double>& xi, DenseMatrix<double>& drdxi) const
444 {
445 throw OomphLibError("dposition() not implemented",
446 OOMPH_CURRENT_FUNCTION,
447 OOMPH_EXCEPTION_LOCATION);
448 }
449 };
450
451} // namespace oomph
452
453#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
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
///////////////////////////////////////////////////////////////////////////// ///////////////////////...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
virtual void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s)=0
Find sub-GeomObject (finite element) and the local coordinate s within it that contains point with gl...
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
A general Finite Element class.
Definition: elements.h:1313
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
void set_nlagrangian_and_ndim(const unsigned &n_lagrangian, const unsigned &n_dim)
Set # of Lagrangian and Eulerian coordinates.
Definition: geom_objects.h:183
unsigned nlagrangian() const
Access function to # of Lagrangian coordinates.
Definition: geom_objects.h:171
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
MeshAsGeomObject()
Empty Constructor.
unsigned Sample_point_container_version
Which version of the sample point container are we using?
unsigned sample_point_container_version() const
Which sample point container is used in locate zeta? (uses enum Sample_Point_Container_Type)
SamplePointContainer * Sample_point_container_pt
Pointer to the sample point container.
void build_it(SamplePointContainerParameters *sample_point_container_parameters_pt)
Helper function to actually build the thing.
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on.
unsigned nelement()
Number of elements in the underlying mesh.
SamplePointContainer * sample_point_container_pt() const
Pointer to the sample point container.
void operator=(const MeshAsGeomObject &)=delete
Broken assignment operator.
OomphCommunicator * Communicator_pt
Communicator.
unsigned ngeom_data() const
How many items of Data does the shape of the object depend on?
Vector< FiniteElement * > Sub_geom_object_pt
Internal storage for the elements that constitute the object.
FiniteElement * finite_element_pt(const unsigned &e)
Return pointer to e-th finite element.
MeshAsGeomObject(Mesh *const &mesh_pt)
Constructor.
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
Find the sub geometric object and local coordinate therein that corresponds to the intrinsic coordina...
MeshAsGeomObject(const MeshAsGeomObject &)=delete
Broken copy constructor.
MeshAsGeomObject(SamplePointContainerParameters *sample_point_container_parameters_pt)
Constructor.
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Parametrised position on object: r(zeta). Evaluated at previous timestep. t=0: current time; t>0: pre...
Mesh * Mesh_pt
Pointer to mesh.
void position(const Vector< double > &zeta, Vector< double > &r) const
Return the position as a function of the intrinsic coordinate zeta. This provides an (expensive!...
void dposition(const Vector< double > &xi, DenseMatrix< double > &drdxi) const
Return the derivative of the position.
Vector< Data * > Geom_data_pt
Vector of pointers to Data items that affects the object's shape.
A general mesh class.
Definition: mesh.h:67
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
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 output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:2027
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed, i.e. if we don't have mpi).
Definition: mesh.h:1600
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
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....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
Mesh * mesh_pt() const
Pointer to mesh from whose FiniteElements sample points are created.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
unsigned Default_sample_point_container_version
Default sample point container type. Must currently be one of UseCGALSamplePointContainer,...
void create_sample_point_container_parameters(Mesh *mesh_pt, SamplePointContainerParameters *&sample_point_container_parameters_pt)
"Factory" for SamplePointContainerParameters of the right type as selected by Default_sample_point_co...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...