// LIC// ==================================================================== // LIC// This file forms part of oomph-lib, the object-oriented, // LIC// multi-physics finite-element library, available // LIC// at http://www.oomph-lib.org. // LIC// // LIC// Copyright (C) 2006-2026 Matthias Heil and Andrew Hazel // LIC// // LIC// This library is free software; you can redistribute it and/or // LIC// modify it under the terms of the GNU Lesser General Public // LIC// License as published by the Free Software Foundation; either // LIC// version 2.1 of the License, or (at your option) any later version. // LIC// // LIC// This library is distributed in the hope that it will be useful, // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU // LIC// Lesser General Public License for more details. // LIC// // LIC// You should have received a copy of the GNU Lesser General Public // LIC// License along with this library; if not, write to the Free Software // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA // LIC// 02110-1301 USA. // LIC// // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk. // LIC// // LIC//==================================================================== // Header file for a relatively simple Quad Meshe #ifndef OOMPH_RECTANGULAR_QUADMESH_HEADER #define OOMPH_RECTANGULAR_QUADMESH_HEADER // Config header #ifdef HAVE_CONFIG_H #include #endif // OOMPH-LIB headers #include "generic/mesh.h" #include "generic/quad_mesh.h" #include "generic/refineable_quad_mesh.h" namespace oomph { //========================================================================== /// RectangularQuadMesh is a two-dimensional mesh of Quad elements with Nx /// elements in the "x" (horizonal) direction and Ny elements in the "y" /// (vertical) direction. Two Constructors are provided. The basic constructor /// assumes that the lower-left-hand corner of the mesh is (0,0) and /// takes only the arguments, Nx, Ny, Xmax and Ymax. The more complex /// constructor takes the additional arguments Xmin and Ymin. /// /// This class is designed to be used as a Base class for more complex /// two dimensional meshes. The virtual functions x_spacing_function() /// and y_spacing_function() may be overloaded to provide arbitrary node /// spacing. The default is uniformly spaced nodes in each direction. /// /// It is also possible to make the solution periodic in the x direction. //=========================================================================== template class RectangularQuadMesh : public virtual QuadMeshBase { protected: // Mesh variables /// Nx: number of elements in x-direction unsigned Nx; /// Ny: number of elements in y-direction unsigned Ny; /// Np: number of (linear) points in the element unsigned Np; /// Minimum value of x coordinate double Xmin; /// Maximum value of x coordinate double Xmax; /// Minimum value of y coordinate double Ymin; /// Maximum value of y coordinate double Ymax; /// Boolean variable used to determine whether the mesh /// is periodic in the x-direction bool Xperiodic; /// Generic mesh construction function: contains all the hard work void build_mesh(TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper); /// Constructor that allows the specification of minimum and maximum /// values of x and y coordinates and does not build the mesh /// This is intend to be used in derived classes that overload the /// spacing functions. THis is scheduled to be changed, however. /// The reason why this MUST be done is because the virtual spacing /// functions cannot be called in the base constructur, because they will /// not have been overloaded yet!! RectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& xmin, const double& xmax, const double& ymin, const double& ymax, const bool& periodic_in_x, const bool& build, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : Nx(nx), Ny(ny), Xmin(xmin), Xmax(xmax), Ymin(ymin), Ymax(ymax), Xperiodic(periodic_in_x) { if (build) { // Call the generic mesh constructor build_mesh(time_stepper_pt); } } public: /// Simple constructor: nx: number of elements in x direction; /// ny: number of elements in y direction; lx, length of domain in x /// direction (0,lx); ly, length of domain in y direction (0,ly) /// Also pass pointer to timestepper (defaults to Steady) RectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : Nx(nx), Ny(ny), Xmin(0.0), Xmax(lx), Ymin(0.0), Ymax(ly), Xperiodic(false) { // Mesh can only be built with 2D Qelements. MeshChecker::assert_geometric_element(2); // Call the generic mesh constructor build_mesh(time_stepper_pt); } /// Constructor that allows the specification of minimum and maximum /// values of x and y coordinates RectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& xmin, const double& xmax, const double& ymin, const double& ymax, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : Nx(nx), Ny(ny), Xmin(xmin), Xmax(xmax), Ymin(ymin), Ymax(ymax), Xperiodic(false) { // Mesh can only be built with 2D Qelements. MeshChecker::assert_geometric_element(2); // Call the generic mesh constructor build_mesh(time_stepper_pt); } /// Simple constructor: nx: number of elements in x direction; /// ny: number of elements in y direction; lx, length of domain in x /// direction (0,lx); ly, length of domain in y direction (0,ly) /// Boolean flag specifies if the mesh is periodic in the x-direction. /// Also pass pointer to timestepper (defaults to Steady) RectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, const bool& periodic_in_x, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : Nx(nx), Ny(ny), Xmin(0.0), Xmax(lx), Ymin(0.0), Ymax(ly), Xperiodic(periodic_in_x) { // Mesh can only be built with 2D Qelements. MeshChecker::assert_geometric_element(2); // Call the generic mesh constructor build_mesh(time_stepper_pt); } /// Constructor that allows the specification of minimum and maximum /// values of x and y coordinates. /// Boolean flag specifies if the mesh is periodic in the x-direction. RectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& xmin, const double& xmax, const double& ymin, const double& ymax, const bool& periodic_in_x, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : Nx(nx), Ny(ny), Xmin(xmin), Xmax(xmax), Ymin(ymin), Ymax(ymax), Xperiodic(periodic_in_x) { // Mesh can only be built with 2D Qelements. MeshChecker::assert_geometric_element(2); // Call the generic mesh constructor build_mesh(time_stepper_pt); } /// Return number of elements in x direction const unsigned& nx() const { // Return the value of Nx return Nx; } /// Return number of elements in y direction const unsigned& ny() const { // Return the value of Ny return Ny; } /// Return the minimum value of x coordinate const double x_min() const { // Return the value of Xmin return Xmin; } /// Return the maximum value of x coordinate const double x_max() const { // Return the value of Xmax return Xmax; } /// Return the minimum value of y coordinate const double y_min() const { // Return the value of Ymin return Ymin; } /// Return the maximum value of y coordinate const double y_max() const { // Return the value of Ymax return Ymax; } /// Reorder the elements: By default they are ordered /// in "horizontal" layers (increasing in x, then in y). This /// function changes this to an ordering in the vertical direction /// (y first, then x). This is more efficient if a frontal solver /// is used and the mesh has more elements in the x than the y direction. /// Can be overloaded in specific derived meshes. virtual void element_reorder(); /// Return the value of the x-coordinate at the node given by the /// local node number (xnode, ynode) in the element (xelement,yelement). /// The description is in a "psudeo" two-dimensional coordinate system, /// so the range of xelement is [0,Nx-1], yelement is [0,Ny-1], and /// that of xnode and ynode is [0,Np-1]. The default is to return /// nodes that are equally spaced in the x coodinate. virtual double x_spacing_function(unsigned xelement, unsigned xnode, unsigned yelement, unsigned ynode) { // Calculate the values of equal increments in nodal values double xstep = (Xmax - Xmin) / ((Np - 1) * Nx); // Return the appropriate value return (Xmin + xstep * ((Np - 1) * xelement + xnode)); } /// Return the value of the y-coordinate at the node given by the /// local node number (xnode, ynode) in the element (xelement,yelement). /// The description is in a "psudeo" two-dimensional coordinate system, /// so the range of xelement is [0,Nx-1], yelement is [0,Ny-1], and /// that of xnode and ynode is [0,Np-1]. The default it to return /// nodes that are equally spaced in the y coordinate. virtual double y_spacing_function(unsigned xelement, unsigned xnode, unsigned yelement, unsigned ynode) { double ystep = (Ymax - Ymin) / ((Np - 1) * Ny); // Return the appropriate value return (Ymin + ystep * ((Np - 1) * yelement + ynode)); } }; //========================================================================== /// Refineable version of the RectangularQuadMesh: A two-dimensional /// mesh of Quad elements with Nx elements in the "x" (horizonal) /// direction and Ny elements in the "y" (vertical) direction. Two /// Constructors are provided. The basic constructor /// assumes that the lower-left-hand corner of the mesh is (0,0) and /// takes only the arguments, Nx, Ny, Xmax and Ymax. The more complex /// constructor takes the additional arguments Xmin and Ymin. /// /// This class is designed to be used as a Base class for more complex /// two dimensional meshes. The virtual functions x_spacing_function() /// and y_spacing_function() may be overloaded to provide arbitrary node /// spacing. The default is uniformly spaced nodes in each direction. //=========================================================================== template class RefineableRectangularQuadMesh : public virtual RectangularQuadMesh, public RefineableQuadMesh { public: /// Simple constructor: nx: number of elements in x direction; /// ny: number of elements in y direction; lx, length of domain in x /// direction (0,lx); ly, length of domain in y direction (0,ly). /// Also pass pointer to timestepper (defaults to Steady) RefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh(nx, ny, lx, ly, time_stepper_pt) { // Nodal positions etc. were created in constructor for // RectangularMesh<...>. Only need to setup quadtree forest this->setup_quadtree_forest(); } /// Simple constructor: nx: number of elements in x direction; /// ny: number of elements in y direction; lx, length of domain in x /// direction (0,lx); ly, length of domain in y direction (0,ly); /// periodic_in_x, periodicity in x. /// Also pass pointer to timestepper (defaults to Steady) RefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, const bool& periodic_in_x, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh( nx, ny, lx, ly, periodic_in_x, time_stepper_pt) { // Nodal positions etc. were created in constructor for // RectangularMesh<...>. Only need to setup quadtree forest this->setup_quadtree_forest(); } /// Constructor that allows the specification of minimum and maximum /// values of x and y coordinates /// Also pass pointer to timestepper (defaults to Steady) RefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& xmin, const double& xmax, const double& ymin, const double& ymax, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh( nx, ny, xmin, xmax, ymin, ymax, time_stepper_pt) { // Nodal positions etc. were created in constructor for // RectangularMesh<...>. Only need to setup quadtree forest this->setup_quadtree_forest(); } /// Constructor that allows the specification of minimum and maximum /// values of x and y coordinates and periodicity /// Also pass pointer to timestepper (defaults to Steady) RefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& xmin, const double& xmax, const double& ymin, const double& ymax, const bool& periodic_in_x, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh( nx, ny, xmin, xmax, ymin, ymax, periodic_in_x, time_stepper_pt) { // Nodal positions etc. were created in constructor for // RectangularMesh<...>. Only need to setup quadtree forest this->setup_quadtree_forest(); } }; ////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// //================================================================ /// Elastic quad mesh with functionality to /// attach traction elements to the specified boundaries. We "upgrade" /// the RectangularQuadMesh to become an /// SolidMesh and equate the Eulerian and Lagrangian coordinates, /// thus making the domain represented by the mesh the stress-free /// configuration. //================================================================ template class ElasticRectangularQuadMesh : public virtual RectangularQuadMesh, public virtual SolidMesh { public: /// Constructor: Build mesh and copy Eulerian coords to Lagrangian /// ones so that the initial configuration is the stress-free one and /// assign boundary coordinates. Origin specifies /// an additional rigid-body displacement. ElasticRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, const Vector& origin, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh(nx, ny, lx, ly, time_stepper_pt) { // Translate the nodes unsigned nnod = nnode(); for (unsigned j = 0; j < nnod; j++) { node_pt(j)->x(0) += origin[0]; node_pt(j)->x(1) += origin[1]; } /// Make the current configuration the undeformed one by /// setting the nodal Lagrangian coordinates to their current /// Eulerian ones set_lagrangian_nodal_coordinates(); // Setup boundary coordinates set_boundary_coordinates(origin); } /// Constructor: Build mesh and copy Eulerian coords to Lagrangian /// ones so that the initial configuration is the stress-free one and /// assign boundary coordinates ElasticRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh(nx, ny, lx, ly, time_stepper_pt) { // No shift Vector origin(2, 0.0); /// Make the current configuration the undeformed one by /// setting the nodal Lagrangian coordinates to their current /// Eulerian ones set_lagrangian_nodal_coordinates(); // Setup boundary coordinates set_boundary_coordinates(origin); } /// Constructor: Build mesh and copy Eulerian coords to Lagrangian /// ones so that the initial configuration is the stress-free one and /// assign boundary coordinates. This includes a boolean flag to specify /// if the mesh is periodic in the x-direction ElasticRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, const bool& periodic_in_x, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh( nx, ny, lx, ly, periodic_in_x, time_stepper_pt) { // No shift Vector origin(2, 0.0); /// Make the current configuration the undeformed one by /// setting the nodal Lagrangian coordinates to their current /// Eulerian ones set_lagrangian_nodal_coordinates(); // Setup boundary coordinates set_boundary_coordinates(origin); } private: /// Setup the boundary coordinates. Vector /// origin specifies the coordinates of the lower left corner of /// the mesh. void set_boundary_coordinates(const Vector& origin) { // 1D vector fo boundary coordinate Vector zeta(1); // Loop over boundaries 0 and 2 where xi_0 is the varying // boundary coordinate for (unsigned b = 0; b < 3; b += 2) { // Number of nodes on those boundaries unsigned n_nod = nboundary_node(b); // Loop over the nodes for (unsigned i = 0; i < n_nod; i++) { // Boundary coordinate varies between 0 and L zeta[0] = boundary_node_pt(b, i)->xi(0) - origin[0]; boundary_node_pt(b, i)->set_coordinates_on_boundary(b, zeta); } set_boundary_coordinate_exists(b); } // Loop over boundaries 1 and 3 where xi_1 is the varying // boundary coordinate for (unsigned b = 1; b < 4; b += 2) { // Number of nodes on those boundaries unsigned n_nod = nboundary_node(b); // Loop over the nodes for (unsigned i = 0; i < n_nod; i++) { // Boundary coordinate varies between +/- H/2 zeta[0] = boundary_node_pt(b, i)->xi(1) - origin[1] - 0.5 * (this->Ymax - this->Ymin); boundary_node_pt(b, i)->set_coordinates_on_boundary(b, zeta); } set_boundary_coordinate_exists(b); } } }; ////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// //================================================================ /// Elastic refineable quad mesh with functionality to /// attach traction elements to the specified boundaries. We "upgrade" /// the RefineableRectangularQuadMesh to become an /// SolidMesh and equate the Eulerian and Lagrangian coordinates, /// thus making the domain represented by the mesh the stress-free /// configuration. We also move the mesh "down" by half the /// the "height" so x=0 is located on the centreline -- appropriate /// for the beam-type problems for which this mesh was developed. //================================================================ template class ElasticRefineableRectangularQuadMesh : public virtual ElasticRectangularQuadMesh, public RefineableQuadMesh { public: /// Constructor: Build mesh and copy Eulerian coords to Lagrangian /// ones so that the initial configuration is the stress-free one and /// assign boundary coordinates (variable Lagrangian coordinates along /// the relevant boundaries). ElasticRefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh(nx, ny, lx, ly, time_stepper_pt), ElasticRectangularQuadMesh(nx, ny, lx, ly, time_stepper_pt) { // Nodal positions etc. were created in base class. // Only need to setup quadtree forest this->setup_quadtree_forest(); } /// Constructor: Build mesh and copy Eulerian coords to Lagrangian /// ones so that the initial configuration is the stress-free one and /// assign boundary coordinates. This includes a boolean flag to specify /// if the mesh is periodic in the x-direction ElasticRefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, const bool& periodic_in_x, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh( nx, ny, lx, ly, periodic_in_x, time_stepper_pt), ElasticRectangularQuadMesh( nx, ny, lx, ly, periodic_in_x, time_stepper_pt) { // Nodal positions etc. were created in base class. // Only need to setup quadtree forest this->setup_quadtree_forest(); } /// Constructor: Build mesh and copy Eulerian coords to Lagrangian /// ones so that the initial configuration is the stress-free one and /// assign boundary coordinates (variable Lagrangian coordinates along /// the relevant boundaries). Origin specifies an additional rigid-body /// displacement. ElasticRefineableRectangularQuadMesh( const unsigned& nx, const unsigned& ny, const double& lx, const double& ly, const Vector& origin, TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper) : RectangularQuadMesh(nx, ny, lx, ly, time_stepper_pt), ElasticRectangularQuadMesh( nx, ny, lx, ly, origin, time_stepper_pt) { // Nodal positions etc. were created in base class. // Only need to setup quadtree forest this->setup_quadtree_forest(); } }; } // namespace oomph #include "rectangular_quadmesh.template.cc" #endif