Public Member Functions | Private Attributes | List of all members
TwoDDGMesh< ELEMENT > Class Template Reference
+ Inheritance diagram for TwoDDGMesh< ELEMENT >:

Public Member Functions

 TwoDDGMesh (const unsigned &Nx, const unsigned &Ny, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
void neighbour_finder (FiniteElement *const &bulk_element_pt, const int &face_index, const Vector< double > &s_bulk, FaceElement *&face_element_pt, Vector< double > &s_face)
 
 TwoDDGMesh (const unsigned &Nx, const unsigned &Ny, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
void neighbour_finder (FiniteElement *const &bulk_element_pt, const int &face_index, const Vector< double > &s_bulk, FaceElement *&face_element_pt, Vector< double > &s_face)
 
 TwoDDGMesh (const unsigned &n_x, const unsigned &n_y, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
void neighbour_finder (FiniteElement *const &bulk_element_pt, const int &face_index, const Vector< double > &s_bulk, FaceElement *&face_element_pt, Vector< double > &s_face)
 
 TwoDDGMesh (const unsigned &Nx, const unsigned &Ny, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
void neighbour_finder (FiniteElement *const &bulk_element_pt, const int &face_index, const Vector< double > &s_bulk, FaceElement *&face_element_pt, Vector< double > &s_face)
 
 TwoDDGMesh (const unsigned &Nx, const unsigned &Ny, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
 
void neighbour_finder (FiniteElement *const &bulk_element_pt, const int &face_index, const Vector< double > &s_bulk, FaceElement *&face_element_pt, Vector< double > &s_face)
 
- Public Member Functions inherited from oomph::DGMesh
 DGMesh ()
 
virtual ~DGMesh ()
 
void setup_face_neighbour_info (const bool &add_face_data_as_external=false)
 
void limit_slopes (SlopeLimiter *const &slope_limiter_pt)
 
- Public Member Functions inherited from oomph::Mesh
 Mesh ()
 Default constructor. More...
 
 Mesh (const Vector< Mesh * > &sub_mesh_pt)
 
void merge_meshes (const Vector< Mesh * > &sub_mesh_pt)
 
virtual void setup_boundary_element_info ()
 
virtual void setup_boundary_element_info (std::ostream &outfile)
 
virtual void reset_boundary_element_info (Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
 Virtual function to perform the reset boundary elements info rutines. More...
 
template<class BULK_ELEMENT >
void doc_boundary_coordinates (const unsigned &b, std::ofstream &the_file)
 
virtual void scale_mesh (const double &factor)
 
 Mesh (const Mesh &dummy)=delete
 Broken copy constructor. More...
 
void operator= (const Mesh &)=delete
 Broken assignment operator. More...
 
virtual ~Mesh ()
 Virtual Destructor to clean up all memory. More...
 
void flush_element_and_node_storage ()
 
void flush_element_storage ()
 
void flush_node_storage ()
 
Node *& node_pt (const unsigned long &n)
 Return pointer to global node n. More...
 
Node * node_pt (const unsigned long &n) const
 Return pointer to global node n (const version) More...
 
GeneralisedElement *& element_pt (const unsigned long &e)
 Return pointer to element e. More...
 
GeneralisedElement * element_pt (const unsigned long &e) const
 Return pointer to element e (const version) More...
 
const Vector< GeneralisedElement * > & element_pt () const
 Return reference to the Vector of elements. More...
 
Vector< GeneralisedElement * > & element_pt ()
 Return reference to the Vector of elements. More...
 
FiniteElement * finite_element_pt (const unsigned &e) const
 
Node *& boundary_node_pt (const unsigned &b, const unsigned &n)
 Return pointer to node n on boundary b. More...
 
Node * boundary_node_pt (const unsigned &b, const unsigned &n) const
 Return pointer to node n on boundary b. More...
 
void set_nboundary (const unsigned &nbound)
 Set the number of boundaries in the mesh. More...
 
void remove_boundary_nodes ()
 Clear all pointers to boundary nodes. More...
 
void remove_boundary_nodes (const unsigned &b)
 
void remove_boundary_node (const unsigned &b, Node *const &node_pt)
 Remove a node from the boundary b. More...
 
void add_boundary_node (const unsigned &b, Node *const &node_pt)
 Add a (pointer to) a node to the b-th boundary. More...
 
void copy_boundary_node_data_from_nodes ()
 
bool boundary_coordinate_exists (const unsigned &i) const
 Indicate whether the i-th boundary has an intrinsic coordinate. More...
 
unsigned long nelement () const
 Return number of elements in the mesh. More...
 
unsigned long nnode () const
 Return number of nodes in the mesh. More...
 
unsigned ndof_types () const
 Return number of dof types in mesh. More...
 
unsigned elemental_dimension () const
 Return number of elemental dimension in mesh. More...
 
unsigned nodal_dimension () const
 Return number of nodal dimension in mesh. More...
 
void add_node_pt (Node *const &node_pt)
 Add a (pointer to a) node to the mesh. More...
 
void add_element_pt (GeneralisedElement *const &element_pt)
 Add a (pointer to) an element to the mesh. More...
 
virtual void node_update (const bool &update_all_solid_nodes=false)
 
virtual void reorder_nodes (const bool &use_old_ordering=true)
 
virtual void get_node_reordering (Vector< Node * > &reordering, const bool &use_old_ordering=true) const
 
template<class BULK_ELEMENT , template< class > class FACE_ELEMENT>
void build_face_mesh (const unsigned &b, Mesh *const &face_mesh_pt)
 
unsigned self_test ()
 Self-test: Check elements and nodes. Return 0 for OK. More...
 
void max_and_min_element_size (double &max_size, double &min_size)
 
double total_size ()
 
void check_inverted_elements (bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
 
void check_inverted_elements (bool &mesh_has_inverted_elements)
 
unsigned check_for_repeated_nodes (const double &epsilon=1.0e-12)
 
Vector< Node * > prune_dead_nodes ()
 
unsigned nboundary () const
 Return number of boundaries. More...
 
unsigned long nboundary_node (const unsigned &ibound) const
 Return number of nodes on a particular boundary. More...
 
FiniteElement * boundary_element_pt (const unsigned &b, const unsigned &e) const
 Return pointer to e-th finite element on boundary b. More...
 
Node * get_some_non_boundary_node () const
 
unsigned nboundary_element (const unsigned &b) const
 Return number of finite elements that are adjacent to boundary b. More...
 
int face_index_at_boundary (const unsigned &b, const unsigned &e) const
 
virtual void dump (std::ofstream &dump_file, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
void dump (const std::string &dump_file_name, const bool &use_old_ordering=true) const
 Dump the data in the mesh into a file for restart. More...
 
virtual void read (std::ifstream &restart_file)
 Read solution from restart file. More...
 
void output_paraview (std::ofstream &file_out, const unsigned &nplot) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt) const
 
void output_fct_paraview (std::ofstream &file_out, const unsigned &nplot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt) const
 
void output (std::ostream &outfile)
 Output for all elements. More...
 
void output (std::ostream &outfile, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output (FILE *file_pt)
 Output for all elements (C-style output) More...
 
void output (FILE *file_pt, const unsigned &nplot)
 Output at f(n_plot) points in each element (C-style output) More...
 
void output (const std::string &output_filename)
 Output for all elements. More...
 
void output (const std::string &output_filename, const unsigned &n_plot)
 Output at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt)
 Output a given Vector function at f(n_plot) points in each element. More...
 
void output_fct (std::ostream &outfile, const unsigned &n_plot, const double &time, FiniteElement::UnsteadyExactSolutionFctPt)
 
void output_boundaries (std::ostream &outfile)
 Output the nodes on the boundaries (into separate tecplot zones) More...
 
void output_boundaries (const std::string &output_filename)
 
void assign_initial_values_impulsive ()
 Assign initial values for an impulsive start. More...
 
void shift_time_values ()
 
void calculate_predictions ()
 
void set_nodal_and_elemental_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void set_mesh_level_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
void set_consistent_pinned_values_for_continuation (ContinuationStorageScheme *const &continuation_stepper_pt)
 Set consistent values for pinned data in continuation. More...
 
bool does_pointer_correspond_to_mesh_data (double *const &parameter_pt)
 Does the double pointer correspond to any mesh data. More...
 
void set_nodal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 Set the timestepper associated with the nodal data in the mesh. More...
 
void set_elemental_internal_time_stepper (TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
 
virtual void compute_norm (double &norm)
 
virtual void compute_norm (Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
 
virtual void compute_error (FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
 Returns the norm of the error and that of the exact solution. More...
 
virtual void compute_error (FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, Vector< double > &error, Vector< double > &norm)
 
bool is_mesh_distributed () const
 Boolean to indicate if Mesh has been distributed. More...
 
OomphCommunicator * communicator_pt () const
 
void delete_all_external_storage ()
 Wipe the storage for all externally-based elements. More...
 

Private Attributes

std::map< std::pair< FiniteElement *, unsigned >, FiniteElement * > Neighbour_map
 
std::map< std::pair< FiniteElement *, int >, FiniteElement * > Neighbour_map
 
std::map< std::pair< FiniteElement *, int >, FaceElement * > Face_element_pt
 

Additional Inherited Members

- Public Types inherited from oomph::Mesh
typedef void(FiniteElement::* SteadyExactSolutionFctPt) (const Vector< double > &x, Vector< double > &soln)
 
typedef void(FiniteElement::* UnsteadyExactSolutionFctPt) (const double &time, const Vector< double > &x, Vector< double > &soln)
 
- Static Public Attributes inherited from oomph::DGMesh
static double FaceTolerance = 1.0e-10
 
- Static Public Attributes inherited from oomph::Mesh
static Steady< 0 > Default_TimeStepper
 The Steady Timestepper. More...
 
static bool Suppress_warning_about_empty_mesh_level_time_stepper_function
 Static boolean flag to control warning about mesh level timesteppers. More...
 
- Protected Member Functions inherited from oomph::Mesh
unsigned long assign_global_eqn_numbers (Vector< double * > &Dof_pt)
 Assign (global) equation numbers to the nodes. More...
 
void describe_dofs (std::ostream &out, const std::string &current_string) const
 
void describe_local_dofs (std::ostream &out, const std::string &current_string) const
 
void assign_local_eqn_numbers (const bool &store_local_dof_pt)
 Assign local equation numbers in all elements. More...
 
void convert_to_boundary_node (Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
 
void convert_to_boundary_node (Node *&node_pt)
 
- Protected Attributes inherited from oomph::Mesh
Vector< Vector< Node * > > Boundary_node_pt
 
bool Lookup_for_elements_next_boundary_is_setup
 
Vector< Vector< FiniteElement * > > Boundary_element_pt
 
Vector< Vector< int > > Face_index_at_boundary
 
Vector< Node * > Node_pt
 Vector of pointers to nodes. More...
 
Vector< GeneralisedElement * > Element_pt
 Vector of pointers to generalised elements. More...
 
std::vector< bool > Boundary_coordinate_exists
 

Constructor & Destructor Documentation

◆ TwoDDGMesh() [1/5]

template<class ELEMENT >
TwoDDGMesh< ELEMENT >::TwoDDGMesh ( const unsigned &  Nx,
const unsigned &  Ny,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline
74  {
75  Vector<double> s_fraction;
76  //Lengths of the mesh
77  double Lx = 1.0;
78  double Ly = 1.0;
79 
80  double llx = -0.5;
81  double lly = -0.5;
82  //Work out the length of each element in the x and y directions
83  //(Assuming uniform spacing)
84  double el_length[2] = {Lx/(double)Nx, Ly/(double)Ny};
85 
86  //loop over the elements in x
87  for(unsigned ex=0;ex<Nx;ex++)
88  {
89  //loop over the element in y
90  for(unsigned ey=0;ey<Ny;ey++)
91  {
92  //Create a new DG element
93  ELEMENT* local_element_pt = new ELEMENT;
94  //Construct nodes and faces
95  local_element_pt->construct_nodes_and_faces(this,time_stepper_pt);
96  //Find the number of nodes
97  unsigned n_node = local_element_pt->nnode();
98  //Find the lower left corner of the element
99  double ll_corner[2] = {llx + ex*el_length[0],lly + ey*el_length[1]};
100  //For each of the nodes set the position
101  for(unsigned n=0;n<n_node;n++)
102  {
103  //Get pointer to the node
104  Node* nod_pt = local_element_pt->node_pt(n);
105  //Get the relative position in local coordinates
106  local_element_pt->local_fraction_of_node(n,s_fraction);
107  //Loop over the coordinates and set the position
108  for(unsigned i=0;i<2;i++)
109  {
110  nod_pt->x(i) = ll_corner[i] + el_length[i]*s_fraction[i];
111  }
112  //Now set the value
113  nod_pt->set_value(0,1.0);
114 
115  //Add each node to the node list
116  Node_pt.push_back(nod_pt);
117  }
118  //Now add the element to the list
119  Element_pt.push_back(local_element_pt);
120  }
121  }
122 
123  //Now loop over all the elements and set up the neighbour map
124  for(unsigned ex=0;ex<Nx;ex++)
125  {
126  for(unsigned ey=0;ey<Ny;ey++)
127  {
128  //Get pointer to the element
129  unsigned element_index = ex*Ny + ey;
130  FiniteElement* local_el_pt = finite_element_pt(element_index);
131 
132  //Storage for indices of neighbours
133  unsigned index[4];
134 
135  //North neighbour (just one up)
136  if(ey < (Ny-1)) {index[0] = element_index + 1;}
137  //If at top, return index at bottom (periodic conditions)
138  else {index[0] = ex*Ny;}
139 
140  //East neighbour (just one across)
141  if(ex < (Nx-1)) {index[1] = element_index + Ny;}
142  //If at side return index at other side (periodic conditions)
143  else {index[1] = ey;}
144 
145  //South neighbour
146  if(ey > 0) {index[2] = element_index - 1;}
147  else {index[2] = element_index + Ny-1;}
148 
149  //West neighbour
150  if(ex > 0) {index[3] = element_index - Ny;}
151  else {index[3] = element_index + (Nx-1)*Ny;}
152 
153  //Now store the details in the mape
154  for(unsigned i=0;i<4;i++)
155  {
156  Neighbour_map[std::make_pair(local_el_pt,i)] =
157  finite_element_pt(index[i]);
158  }
159  }
160  }
161 
162  //Setup the connections between the neighbouring faces
163  this->setup_face_neighbour_info();
164 
165  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
TwoDDGMesh::Neighbour_map
std::map< std::pair< FiniteElement *, unsigned >, FiniteElement * > Neighbour_map
Definition: two_d_advection.cc:68
oomph::DGMesh::setup_face_neighbour_info
void setup_face_neighbour_info(const bool &add_face_data_as_external=false)
Definition: dg_elements.h:489
oomph::Data::set_value
void set_value(const unsigned &i, const double &value_)
Definition: nodes.h:271
oomph::FiniteElement
Definition: elements.h:1313
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Definition: mesh.h:473
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
oomph::Node
Definition: nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
GlobalParameters::Nx
unsigned Nx
Number of elements in each direction (used by SimpleCubicMesh)
Definition: structured_cubic_point_source.cc:114
GlobalParameters::Ny
unsigned Ny
Definition: structured_cubic_point_source.cc:115
Global_Parameters::Ly
double Ly
Length of domain in y direction.
Definition: periodic_load.cc:58
Global_Parameters::Lx
double Lx
Length of domain in x direction.
Definition: periodic_load.cc:55

References i, Global_Parameters::Lx, Global_Parameters::Ly, n, GlobalParameters::Nx, GlobalParameters::Ny, oomph::Data::set_value(), and oomph::Node::x().

◆ TwoDDGMesh() [2/5]

template<class ELEMENT >
TwoDDGMesh< ELEMENT >::TwoDDGMesh ( const unsigned &  Nx,
const unsigned &  Ny,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline
78  {
79  Vector<double> s_fraction;
80  //Lengths of the mesh
81  double R_min = 1.0;
82  double R_max = 4.0;
83 
84  const double pi = MathematicalConstants::Pi;
85 
86  //Work out the length of each element in the x and y directions
87  //(Assuming uniform spacing)
88  double el_length[2] = {(R_max-R_min)/(double)Nx, (2.0*pi)/(double)Ny};
89 
90  //Vector of booleans for boundary information
91  std::vector<bool> boundary_info;
92 
93  //loop over the elements in x
94  for(unsigned ex=0;ex<Nx;ex++)
95  {
96  //loop over the element in y
97  for(unsigned ey=0;ey<Ny;ey++)
98  {
99  //Create a new DG element
100  ELEMENT* local_element_pt = new ELEMENT;
101  //Find the number of nodes
102  const unsigned n_node = local_element_pt->nnode();
103  //Have we constructed
104  bool constructed = false;
105 
106  //If we are on a boundary
107  if(ex==0)
108  {
109  const unsigned n_p = local_element_pt->nnode_1d();
110  boundary_info.resize(n_node);
111  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
112  //The left hand side is on a boundary
113  for(unsigned n=0;n<n_p;n++)
114  {
115  boundary_info[n*n_p] = true;
116  }
117 
118  //Lower left
119  if(ey==0)
120  {
121  for(unsigned n=1;n<n_p;n++) {boundary_info[n] = true;}
122  }
123 
124  //Upper left
125  if(ey==Ny-1)
126  {
127  for(unsigned n=1;n<n_p;n++) {boundary_info[n_p*(n_p-1) + n] = true;}
128  }
129 
130  //Now construct
131  local_element_pt->construct_boundary_nodes_and_faces(
132  this, boundary_info, time_stepper_pt);
133 
134  constructed = true;
135  }
136  //Right boundary
137  else if(ex==Nx-1)
138  {
139  const unsigned n_p = local_element_pt->nnode_1d();
140  boundary_info.resize(n_node);
141  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
142  //The right-hand side is on a boundary
143  for(unsigned n=0;n<n_p;n++)
144  {
145  boundary_info[n_p-1 + n*n_p] = true;
146  }
147 
148  //Lower right
149  if(ey==0)
150  {
151  for(unsigned n=0;n<n_p-1;n++) {boundary_info[n] = true;}
152  }
153 
154  //Upper right
155  if(ey==Ny-1)
156  {
157  for(unsigned n=0;n<n_p-1;n++)
158  {boundary_info[n_p*(n_p-1) + n] = true;}
159  }
160 
161  //Now construct
162  local_element_pt->construct_boundary_nodes_and_faces(
163  this, boundary_info, time_stepper_pt);
164 
165  constructed = true;
166  }
167  //Otherwise it's in the middle
168  else
169  {
170  //If on the bottom or top
171  if((ey==0) || (ey==Ny-1))
172  {
173  const unsigned n_p = local_element_pt->nnode_1d();
174  boundary_info.resize(n_node);
175  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
176 
177  //If the bottom
178  if(ey==0)
179  {
180  for(unsigned n=0;n<n_p;n++)
181  {
182  boundary_info[n] = true;
183  }
184  }
185 
186  //If on the top
187  if(ey==Ny-1)
188  {
189  for(unsigned n=0;n<n_p;n++)
190  {boundary_info[n_p*(n_p-1) + n] = true;}
191  }
192 
193 
194  //Now construct
195  local_element_pt->construct_boundary_nodes_and_faces(
196  this, boundary_info, time_stepper_pt);
197 
198  constructed = true;
199  }
200  }
201 
202  //If we haven't already constructed
203  if(!constructed)
204  {
205  //Construct nodes and faces
206  local_element_pt->construct_nodes_and_faces(this,time_stepper_pt);
207  }
208 
209  //Find the lower left corner of the element
210  double ll_corner[2] = {R_min + ex*el_length[0],ey*el_length[1]};
211  //For each of the nodes set the position
212  for(unsigned n=0;n<n_node;n++)
213  {
214  //Get pointer to the node
215  Node* nod_pt = local_element_pt->node_pt(n);
216  //Get the relative position in local coordinates
217  local_element_pt->local_fraction_of_node(n,s_fraction);
218  //Loop over the coordinates and set the position
219  double r = ll_corner[0] + el_length[0]*s_fraction[0];
220  double theta = ll_corner[1] + el_length[1]*s_fraction[1];
221  nod_pt->x(0) = r*cos(theta);
222  nod_pt->x(1) = r*sin(theta);
223 
224  //Add each node to the node list
225  Node_pt.push_back(nod_pt);
226  }
227  //Now add the element to the list
228  Element_pt.push_back(local_element_pt);
229  }
230  }
231 
232  //Set up the boundary nodes
233  this->set_nboundary(2);
234 
235  //Now let's add the boundary nodes
236  //Left-hand and right-hand sides
237  for(unsigned e=0;e<Ny;e++)
238  {
239  //Left
240  FiniteElement* elem_pt = this->finite_element_pt(e);
241  unsigned n_p = elem_pt->nnode_1d();
242  for(unsigned n=0;n<n_p;n++)
243  {
244  this->add_boundary_node(0,elem_pt->node_pt(n_p*n));
245  }
246  //Right
247  elem_pt = this->finite_element_pt(Ny*(Nx-1)+e);
248  n_p = elem_pt->nnode_1d();
249  for(unsigned n=0;n<n_p;n++)
250  {
251  this->add_boundary_node(1,elem_pt->node_pt(n_p-1 + n_p*n));
252  }
253  }
254 
255  //Now loop over all the elements and set up the neighbour map
256  for(unsigned ex=0;ex<Nx;ex++)
257  {
258  for(unsigned ey=0;ey<Ny;ey++)
259  {
260  //Get pointer to the element
261  unsigned element_index = ex*Ny + ey;
262  FiniteElement* local_el_pt = finite_element_pt(element_index);
263 
264  //Storage for indices of neighbours
265  unsigned index[4];
266 
267  //North neighbour (just one up)
268  if(ey < (Ny-1)) {index[0] = element_index + 1;}
269  //If at top, return index at bottom (periodic conditions)
270  else {index[0] = ex*Ny;}
271 
272  //East neighbour (just one across)
273  if(ex < (Nx-1)) {index[1] = element_index + Ny;}
274  //If at side return index at other side (periodic conditions)
275  else {index[1] = ey;}
276 
277  //South neighbour
278  if(ey > 0) {index[2] = element_index - 1;}
279  else {index[2] = element_index + Ny-1;}
280 
281  //West neighbour
282  if(ex > 0) {index[3] = element_index - Ny;}
283  else {index[3] = element_index + (Nx-1)*Ny;}
284 
285  //Now store the details in the mape
286  for(unsigned i=0;i<4;i++)
287  {
288  Neighbour_map[std::make_pair(local_el_pt,i)] =
289  finite_element_pt(index[i]);
290  }
291  }
292  }
293 
294  //Now set up the neighbour information
295  /*const unsigned n_element = this->nelement();
296  for(unsigned e=0;e<n_element;e++)
297  {
298  dynamic_cast<ELEMENT*>(this->element_pt(e))->setup_face_neighbour_info();
299  }*/
300 
301  this->setup_face_neighbour_info();
302  }
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
sin
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::nnode_1d
virtual unsigned nnode_1d() const
Definition: elements.h:2218
oomph::Mesh::add_boundary_node
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
oomph::Mesh::set_nboundary
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
BiharmonicTestFunctions2::Pi
double Pi
Definition: two_d_biharmonic.cc:235
BiharmonicTestFunctions2::theta
double theta
Definition: two_d_biharmonic.cc:236
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
constants::pi
const Mdouble pi
Definition: ExtendedMath.h:23

References cos(), e(), i, n, oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), GlobalParameters::Nx, GlobalParameters::Ny, constants::pi, oomph::MathematicalConstants::Pi, UniformPSDSelfTest::r, sin(), BiharmonicTestFunctions2::theta, and oomph::Node::x().

◆ TwoDDGMesh() [3/5]

template<class ELEMENT >
TwoDDGMesh< ELEMENT >::TwoDDGMesh ( const unsigned &  n_x,
const unsigned &  n_y,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline
175  {
176  Vector<double> s_fraction;
177  //Lengths of the mesh
178  double Lx = 2.4;
179  double Ly = 0.8;
180 
181  unsigned Nx = 4*n_x, Ny = 4*n_y;
182 
183  //Work out the length of each element in the x and y directions
184  //(Assuming uniform spacing)
185  //Create the flow bit
186  double el_length[2] = {Lx/(double)Nx, Ly/(double)Ny};
187 
188  //Vector of booleans for boundary information
189  std::vector<bool> boundary_info;
190 
191  unsigned Nx1 = n_x;
192  unsigned Ny1 = 5*n_y;
193  double llx1 = 0.0;
194  double lly1 = 0.0;
195 
196  //loop over the elements in x
197  for(unsigned ex=0;ex<Nx1;ex++)
198  {
199  //loop over the element in y
200  for(unsigned ey=0;ey<Ny1;ey++)
201  {
202  //Create a new DG element
203  ELEMENT* local_element_pt = new ELEMENT;
204  //Find the number of nodes
205  const unsigned n_node = local_element_pt->nnode();
206  //Have we constructed
207  bool constructed = false;
208  //Left boundary
209  if(ex==0)
210  {
211  const unsigned n_p = local_element_pt->nnode_1d();
212  boundary_info.resize(n_node);
213  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
214  //The left hand side is on a boundary
215  for(unsigned n=0;n<n_p;n++)
216  {
217  boundary_info[n*n_p] = true;
218  }
219 
220  //Lower left
221  if(ey==0)
222  {
223  for(unsigned n=1;n<n_p;n++) {boundary_info[n] = true;}
224  }
225 
226  //Upper left
227  if(ey==Ny1-1)
228  {
229  for(unsigned n=1;n<n_p;n++) {boundary_info[n_p*(n_p-1) + n] = true;}
230  }
231 
232  //Now construct
233  local_element_pt->construct_boundary_nodes_and_faces(
234  this, boundary_info, time_stepper_pt);
235 
236  constructed = true;
237  }
238  //Right boundary
239  else if(ex==Nx1-1)
240  {
241  const unsigned n_p = local_element_pt->nnode_1d();
242  boundary_info.resize(n_node);
243  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
244 
245  //Bottom right
246  if(ey==0)
247  {
248  //The bottom is a boundary
249  for(unsigned n=0;n<n_p-1;n++) {boundary_info[n] = true;}
250  }
251 
252 
253  //The right-hand side is on a boundary in the early bit
254  //of the step
255  if(ey < n_y)
256  {
257  for(unsigned n=0;n<n_p;n++)
258  {
259  boundary_info[n_p-1 + n*n_p] = true;
260  }
261  }
262 
263  //If top right
264  if(ey==Ny1-1)
265  {
266  //The top is a boundary
267  for(unsigned n=0;n<n_p;n++)
268  {boundary_info[n_p*(n_p-1) + n] = true;}
269  }
270 
271  local_element_pt->construct_boundary_nodes_and_faces(
272  this, boundary_info, time_stepper_pt);
273 
274  constructed = true;
275  }
276  //Otherwise it's in the middle
277  else
278  {
279  //If on the bottom or top
280  if((ey==0) || (ey==Ny1-1))
281  {
282  const unsigned n_p = local_element_pt->nnode_1d();
283  boundary_info.resize(n_node);
284  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
285 
286  //If the bottom
287  if(ey==0)
288  {
289  for(unsigned n=0;n<n_p;n++)
290  {
291  boundary_info[n] = true;
292  }
293  }
294 
295  //If on the top
296  if(ey==Ny1-1)
297  {
298  for(unsigned n=0;n<n_p;n++)
299  {boundary_info[n_p*(n_p-1) + n] = true;}
300  }
301 
302 
303  //Now construct
304  local_element_pt->construct_boundary_nodes_and_faces(
305  this, boundary_info, time_stepper_pt);
306 
307  constructed = true;
308  }
309  }
310 
311  //If we haven't already constructed
312  if(!constructed)
313  {
314  //Construct nodes and faces
315  local_element_pt->construct_nodes_and_faces(this,time_stepper_pt);
316  }
317 
318  //Find the lower left corner of the element
319  double ll_corner[2] = {llx1 + ex*el_length[0],lly1 + ey*el_length[1]};
320  //For each of the nodes set the position
321  for(unsigned n=0;n<n_node;n++)
322  {
323  //Get pointer to the node
324  Node* nod_pt = local_element_pt->node_pt(n);
325  //Get the relative position in local coordinates
326  local_element_pt->local_fraction_of_node(n,s_fraction);
327  //Loop over the coordinates and set the position
328  for(unsigned i=0;i<2;i++)
329  {
330  nod_pt->x(i) = ll_corner[i] + el_length[i]*s_fraction[i];
331  }
332 
333  //Add each node to the node list
334  Node_pt.push_back(nod_pt);
335  }
336  //Now add the element to the list
337  Element_pt.push_back(local_element_pt);
338  }
339  }
340 
341 
342  //Now do the main bit of the mesh
343  double llx = 0.6;
344  double lly = 0.2;
345 
346  //loop over the elements in x
347  for(unsigned ex=0;ex<Nx;ex++)
348  {
349  //loop over the element in y
350  for(unsigned ey=0;ey<Ny;ey++)
351  {
352  //Create a new DG element
353  ELEMENT* local_element_pt = new ELEMENT;
354  //Find the number of nodes
355  const unsigned n_node = local_element_pt->nnode();
356  //Have we constructed
357  bool constructed = false;
358  //Right boundary
359  if(ex==Nx-1)
360  {
361  const unsigned n_p = local_element_pt->nnode_1d();
362  boundary_info.resize(n_node);
363  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
364  //The right-hand side is on a boundary
365  for(unsigned n=0;n<n_p;n++)
366  {
367  boundary_info[n_p-1 + n*n_p] = true;
368  }
369 
370  //Lower right
371  if(ey==0)
372  {
373  for(unsigned n=0;n<n_p-1;n++) {boundary_info[n] = true;}
374  }
375 
376  //Upper right
377  if(ey==Ny-1)
378  {
379  for(unsigned n=0;n<n_p-1;n++)
380  {boundary_info[n_p*(n_p-1) + n] = true;}
381  }
382 
383  //Now construct
384  local_element_pt->construct_boundary_nodes_and_faces(
385  this, boundary_info, time_stepper_pt);
386 
387  constructed = true;
388  }
389  //Otherwise it's in the middle
390  else
391  {
392  //If on the bottom or top
393  if((ey==0) || (ey==Ny-1))
394  {
395  const unsigned n_p = local_element_pt->nnode_1d();
396  boundary_info.resize(n_node);
397  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
398 
399  //If the bottom
400  if(ey==0)
401  {
402  for(unsigned n=0;n<n_p;n++)
403  {
404  boundary_info[n] = true;
405  }
406  }
407 
408  //If on the top
409  if(ey==Ny-1)
410  {
411  for(unsigned n=0;n<n_p;n++)
412  {boundary_info[n_p*(n_p-1) + n] = true;}
413  }
414 
415 
416  //Now construct
417  local_element_pt->construct_boundary_nodes_and_faces(
418  this, boundary_info, time_stepper_pt);
419 
420  constructed = true;
421  }
422  }
423 
424  //If we haven't already constructed
425  if(!constructed)
426  {
427  //Construct nodes and faces
428  local_element_pt->construct_nodes_and_faces(this,time_stepper_pt);
429  }
430 
431  //Find the lower left corner of the element
432  double ll_corner[2] = {llx + ex*el_length[0],lly + ey*el_length[1]};
433  //For each of the nodes set the position
434  for(unsigned n=0;n<n_node;n++)
435  {
436  //Get pointer to the node
437  Node* nod_pt = local_element_pt->node_pt(n);
438  //Get the relative position in local coordinates
439  local_element_pt->local_fraction_of_node(n,s_fraction);
440  //Loop over the coordinates and set the position
441  for(unsigned i=0;i<2;i++)
442  {
443  nod_pt->x(i) = ll_corner[i] + el_length[i]*s_fraction[i];
444  }
445 
446  //Add each node to the node list
447  Node_pt.push_back(nod_pt);
448  }
449  //Now add the element to the list
450  Element_pt.push_back(local_element_pt);
451  }
452  }
453 
454 
455  //Set up the boundary nodes
456  this->set_nboundary(4);
457 
458  //Now let's add the boundary nodes
459  //Left-hand side
460  for(unsigned e=0;e<Ny1;e++)
461  {
462  FiniteElement* const elem_pt = this->finite_element_pt(e);
463  unsigned n_p = elem_pt->nnode_1d();
464  for(unsigned n=0;n<n_p;n++)
465  {
466  this->add_boundary_node(3,elem_pt->node_pt(n_p*n));
467  }
468  }
469 
470  //Right hand side
471  for(unsigned e=0;e<Ny;e++)
472  {
473  FiniteElement* const elem_pt =
474  this->finite_element_pt(Ny1*Nx1 + Ny*(Nx-1)+e);
475  unsigned n_p = elem_pt->nnode_1d();
476  for(unsigned n=0;n<n_p;n++)
477  {
478  this->add_boundary_node(1,elem_pt->node_pt(n_p-1 + n_p*n));
479  }
480  }
481 
482 
483 
484  //Top
485  {
486  //First bit
487  for(unsigned e=0;e<Nx1;e++)
488  {
489  FiniteElement* const elem_pt = this->finite_element_pt(Ny1-1 + Ny1*e);
490  unsigned n_p = elem_pt->nnode_1d();
491  for(unsigned n=0;n<n_p;n++)
492  {
493  this->add_boundary_node(2,elem_pt->node_pt((n_p-1)*n_p + n));
494  }
495 
496  //Now add the reflecting element for the north face
497  Face_element_pt[std::make_pair(elem_pt,0)] =
498  new DGEulerFaceReflectionElement<ELEMENT>(elem_pt,2);
499  }
500 
501  //Second bit
502  for(unsigned e=0;e<Nx;e++)
503  {
504  FiniteElement* const elem_pt =
505  this->finite_element_pt(Nx1*Ny1 + Ny-1 + Ny*e);
506  unsigned n_p = elem_pt->nnode_1d();
507  for(unsigned n=0;n<n_p;n++)
508  {
509  this->add_boundary_node(2,elem_pt->node_pt((n_p-1)*n_p + n));
510  }
511 
512  //Now add the reflecting element for the north face
513  Face_element_pt[std::make_pair(elem_pt,0)] =
514  new DGEulerFaceReflectionElement<ELEMENT>(elem_pt,2);
515  }
516  }
517 
518  //Bottom
519  {
520  //First bit
521  for(unsigned e=0;e<Nx1;e++)
522  {
523  //Bottom
524  FiniteElement* const elem_pt = this->finite_element_pt(Ny1*e);
525  unsigned n_p = elem_pt->nnode_1d();
526  for(unsigned n=0;n<n_p;n++)
527  {
528  this->add_boundary_node(0,elem_pt->node_pt(n));
529  }
530 
531  //Now add the reflecting element for the south face
532  Face_element_pt[std::make_pair(elem_pt,2)] =
533  new DGEulerFaceReflectionElement<ELEMENT>(elem_pt,-2);
534  }
535 
536  //Add the right-hand of the lower elements
537  for(unsigned e=0;e<n_y;e++)
538  {
539  FiniteElement* const elem_pt = this->finite_element_pt(Ny1*(Nx1-1) + e);
540  unsigned n_p = elem_pt->nnode_1d();
541  for(unsigned n=0;n<n_p;n++)
542  {
543  this->add_boundary_node(0,elem_pt->node_pt(n_p-1 + n_p*n));
544  }
545  //Now add the reflecting element for the east face
546  Face_element_pt[std::make_pair(elem_pt,1)] =
547  new DGEulerFaceReflectionElement<ELEMENT>(elem_pt,1);
548 
549  }
550 
551  //Second bit
552  for(unsigned e=0;e<Nx;e++)
553  {
554  //Bottom
555  FiniteElement* const elem_pt =
556  this->finite_element_pt(Nx1*Ny1 + Ny*e);
557  unsigned n_p = elem_pt->nnode_1d();
558  for(unsigned n=0;n<n_p;n++)
559  {
560  this->add_boundary_node(0,elem_pt->node_pt(n));
561  }
562  //Now add the reflecting element for the south face
563  Face_element_pt[std::make_pair(elem_pt,2)] =
564  new DGEulerFaceReflectionElement<ELEMENT>(elem_pt,-2);
565  }
566  }
567 
568  /* for(unsigned b=0;b<4;b++)
569  {
570  std::cout << "Boundary : " << b << "\n";
571  unsigned nbound = this->nboundary_node(b);
572  for(unsigned n=0;n<nbound;n++)
573  {
574  std::cout << "( " << this->boundary_node_pt(b,n)->x(0) << ", "
575  << this->boundary_node_pt(b,n)->x(1) << ")\n";
576  }
577  }*/
578 
579  //How many faces have we constructed
580  std::cout << Face_element_pt.size() << "\n";
581 
582 
583  //Now loop over all the elements and set up the neighbour map
584  for(unsigned ex=0;ex<Nx;ex++)
585  {
586  for(unsigned ey=0;ey<Ny;ey++)
587  {
588  //Get pointer to the element
589  unsigned element_index = Nx1*Ny1 + ex*Ny + ey;
590  FiniteElement* local_el_pt = finite_element_pt(element_index);
591 
592  //Storage for indices of neighbours
593  int index[4];
594 
595  //North neighbour (just one up)
596  if(ey < (Ny-1)) {index[0] = element_index + 1;}
597  //If at top (no neighbour)
598  else {index[0] = -1;}
599 
600  //East neighbour (just one across)
601  if(ex < (Nx-1)) {index[1] = element_index + Ny;}
602  //If at side (no neighbour)
603  else {index[1] = -1;}
604 
605  //South neighbour
606  if(ey > 0) {index[2] = element_index - 1;}
607  else {index[2] = -1;}
608 
609  //West neighbour
610  /*if(ex > 0)*/ {index[3] = element_index - Ny;}
611  //else {index[3] = -1;}
612 
613  //Now store the details in the mape
614  for(unsigned i=0;i<4;i++)
615  {
616  if(index[i] != -1)
617  {
618  Neighbour_map[std::make_pair(local_el_pt,i)] =
619  finite_element_pt(index[i]);
620  }
621  }
622  }
623  }
624 
625  //Now sort out the maps for the first bit
626  for(unsigned ex=0;ex<Nx1;ex++)
627  {
628  for(unsigned ey=0;ey<Ny1;ey++)
629  {
630  //Get pointer to the element
631  unsigned element_index = ex*Ny1 + ey;
632  FiniteElement* local_el_pt = finite_element_pt(element_index);
633 
634  //Storage for indices of neighbours
635  int index[4];
636 
637  //North neighbour (just one up)
638  if(ey < (Ny1-1)) {index[0] = element_index + 1;}
639  //If at top (no neighbour)
640  else {index[0] = -1;}
641 
642  //East neighbour (just one across)
643  if(ex < (Nx1-1)) {index[1] = element_index + Ny1;}
644  //If at side (no neighbour)
645  else
646  {
647  if(ey >= n_y) {index[1] = element_index + Ny;}
648  else {index[1] = -1;}
649  }
650 
651  //South neighbour
652  if(ey > 0) {index[2] = element_index - 1;}
653  else {index[2] = -1;}
654 
655  //West neighbour
656  if(ex > 0) {index[3] = element_index - Ny1;}
657  else {index[3] = -1;}
658 
659  //Now store the details in the mape
660  for(unsigned i=0;i<4;i++)
661  {
662  if(index[i] != -1)
663  {
664  Neighbour_map[std::make_pair(local_el_pt,i)] =
665  finite_element_pt(index[i]);
666  }
667  }
668  }
669  }
670 
671 
672  //Now set up the neighbour information
673  /*const unsigned n_element = this->nelement();
674  for(unsigned e=0;e<n_element;e++)
675  {
676  dynamic_cast<ELEMENT*>(this->element_pt(e))->setup_face_neighbour_info();
677  }*/
678  this->setup_face_neighbour_info();
679 }
TwoDDGMesh::Face_element_pt
std::map< std::pair< FiniteElement *, int >, FaceElement * > Face_element_pt
Definition: ff_step.cc:169
oomph::DGEulerFaceReflectionElement
Definition: euler_elements.h:698

References e(), i, Global_Parameters::Lx, Global_Parameters::Ly, n, oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), GlobalParameters::Nx, GlobalParameters::Ny, and oomph::Node::x().

◆ TwoDDGMesh() [4/5]

template<class ELEMENT >
TwoDDGMesh< ELEMENT >::TwoDDGMesh ( const unsigned &  Nx,
const unsigned &  Ny,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline
94  {
95  Vector<double> s_fraction;
96  //Lengths of the mesh
97  double Lx = 1.0;
98  double Ly = 1.0;
99 
100  double llx = -0.5;
101  double lly = 0.0;
102  //Angle (in degrees of ramp)
103  double theta = 10.0;
104  const double pi = MathematicalConstants::Pi;
105  //Work out the length of each element in the x and y directions
106  //(Assuming uniform spacing)
107  double el_length[2] = {Lx/(double)Nx, Ly/(double)Ny};
108 
109  //Vector of booleans for boundary information
110  std::vector<bool> boundary_info;
111 
112  //loop over the elements in x
113  for(unsigned ex=0;ex<Nx;ex++)
114  {
115  //loop over the element in y
116  for(unsigned ey=0;ey<Ny;ey++)
117  {
118  //Create a new DG element
119  ELEMENT* local_element_pt = new ELEMENT;
120  //Find the number of nodes
121  const unsigned n_node = local_element_pt->nnode();
122  //Have we constructed
123  bool constructed = false;
124 
125  //If we are on a boundary
126  if(ex==0)
127  {
128  const unsigned n_p = local_element_pt->nnode_1d();
129  boundary_info.resize(n_node);
130  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
131  //The left hand side is on a boundary
132  for(unsigned n=0;n<n_p;n++)
133  {
134  boundary_info[n*n_p] = true;
135  }
136 
137  //Lower left
138  if(ey==0)
139  {
140  for(unsigned n=1;n<n_p;n++) {boundary_info[n] = true;}
141  }
142 
143  //Upper left
144  if(ey==Ny-1)
145  {
146  for(unsigned n=1;n<n_p;n++) {boundary_info[n_p*(n_p-1) + n] = true;}
147  }
148 
149  //Now construct
150  local_element_pt->construct_boundary_nodes_and_faces(
151  this, boundary_info, time_stepper_pt);
152 
153  constructed = true;
154  }
155  //Right boundary
156  else if(ex==Nx-1)
157  {
158  const unsigned n_p = local_element_pt->nnode_1d();
159  boundary_info.resize(n_node);
160  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
161  //The right-hand side is on a boundary
162  for(unsigned n=0;n<n_p;n++)
163  {
164  boundary_info[n_p-1 + n*n_p] = true;
165  }
166 
167  //Lower right
168  if(ey==0)
169  {
170  for(unsigned n=0;n<n_p-1;n++) {boundary_info[n] = true;}
171  }
172 
173  //Upper right
174  if(ey==Ny-1)
175  {
176  for(unsigned n=0;n<n_p-1;n++)
177  {boundary_info[n_p*(n_p-1) + n] = true;}
178  }
179 
180  //Now construct
181  local_element_pt->construct_boundary_nodes_and_faces(
182  this, boundary_info, time_stepper_pt);
183 
184  constructed = true;
185  }
186  //Otherwise it's in the middle
187  else
188  {
189  //If on the bottom or top
190  if((ey==0) || (ey==Ny-1))
191  {
192  const unsigned n_p = local_element_pt->nnode_1d();
193  boundary_info.resize(n_node);
194  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
195 
196  //If the bottom
197  if(ey==0)
198  {
199  for(unsigned n=0;n<n_p;n++)
200  {
201  boundary_info[n] = true;
202  }
203  }
204 
205  //If on the top
206  if(ey==Ny-1)
207  {
208  for(unsigned n=0;n<n_p;n++)
209  {boundary_info[n_p*(n_p-1) + n] = true;}
210  }
211 
212 
213  //Now construct
214  local_element_pt->construct_boundary_nodes_and_faces(
215  this, boundary_info, time_stepper_pt);
216 
217  constructed = true;
218  }
219  }
220 
221  //If we haven't already constructed
222  if(!constructed)
223  {
224  //Construct nodes and faces
225  local_element_pt->construct_nodes_and_faces(this,time_stepper_pt);
226  }
227 
228  //Find the lower left corner of the element
229  double ll_corner[2] = {llx + ex*el_length[0],lly + ey*el_length[1]};
230  //For each of the nodes set the position
231  for(unsigned n=0;n<n_node;n++)
232  {
233  //Get pointer to the node
234  Node* nod_pt = local_element_pt->node_pt(n);
235  //Get the relative position in local coordinates
236  local_element_pt->local_fraction_of_node(n,s_fraction);
237  //Loop over the coordinates and set the position
238  for(unsigned i=0;i<2;i++)
239  {
240  nod_pt->x(i) = ll_corner[i] + el_length[i]*s_fraction[i];
241  }
242 
243  //Add the ramp part
244  if(nod_pt->x(0) >= 0.0)
245  {
246  //Work out position of ramp
247  double y_min = nod_pt->x(0)*tan(theta*pi/180.0);
248  //Now rescale the height
249  double frac = nod_pt->x(1)*(Ly-y_min)/Ly;
250  //And add it
251  nod_pt->x(1) = y_min + frac;
252  }
253 
254  //Add each node to the node list
255  Node_pt.push_back(nod_pt);
256  }
257  //Now add the element to the list
258  Element_pt.push_back(local_element_pt);
259  }
260  }
261 
262  //Set up the boundary nodes
263  this->set_nboundary(4);
264 
265  //Now let's add the boundary nodes
266  //Left-hand and right-hand sides
267  for(unsigned e=0;e<Ny;e++)
268  {
269  //Left
270  FiniteElement* elem_pt = this->finite_element_pt(e);
271  unsigned n_p = elem_pt->nnode_1d();
272  for(unsigned n=0;n<n_p;n++)
273  {
274  this->add_boundary_node(3,elem_pt->node_pt(n_p*n));
275  }
276  //Right
277  elem_pt = this->finite_element_pt(Ny*(Nx-1)+e);
278  n_p = elem_pt->nnode_1d();
279  for(unsigned n=0;n<n_p;n++)
280  {
281  this->add_boundary_node(1,elem_pt->node_pt(n_p-1 + n_p*n));
282  }
283  }
284 
285  //Top and bottom
286  for(unsigned e=0;e<Nx;e++)
287  {
288  //Bottom
289  FiniteElement* elem_pt = this->finite_element_pt(Ny*e);
290  unsigned n_p = elem_pt->nnode_1d();
291  for(unsigned n=0;n<n_p;n++)
292  {
293  this->add_boundary_node(0,elem_pt->node_pt(n));
294  }
295 
296  //Now add the reflecting element for the south face
297  Face_element_pt[std::make_pair(elem_pt,2)] =
298  new DGEulerFaceReflectionElement<ELEMENT>(elem_pt,-2);
299 
300  //Top
301  elem_pt = this->finite_element_pt(Ny-1 + Ny*e);
302  n_p = elem_pt->nnode_1d();
303  for(unsigned n=0;n<n_p;n++)
304  {
305  this->add_boundary_node(2,elem_pt->node_pt((n_p-1)*n_p + n));
306  }
307  }
308 
309  //Now loop over all the elements and set up the neighbour map
310  for(unsigned ex=0;ex<Nx;ex++)
311  {
312  for(unsigned ey=0;ey<Ny;ey++)
313  {
314  //Get pointer to the element
315  unsigned element_index = ex*Ny + ey;
316  FiniteElement* local_el_pt = finite_element_pt(element_index);
317 
318  //Storage for indices of neighbours
319  int index[4];
320 
321  //North neighbour (just one up)
322  if(ey < (Ny-1)) {index[0] = element_index + 1;}
323  //If at top, no neighbour
324  else {index[0] = -1;}
325 
326  //East neighbour (just one across)
327  if(ex < (Nx-1)) {index[1] = element_index + Ny;}
328  //If at side no neighbour
329  else {index[1] = -1;}
330 
331  //South neighbour
332  if(ey > 0) {index[2] = element_index - 1;}
333  else {index[2] = -1;}
334 
335  //West neighbour
336  if(ex > 0) {index[3] = element_index - Ny;}
337  else {index[3] = -1;}
338 
339  //Now store the details in the mape
340  for(unsigned i=0;i<4;i++)
341  {
342  if(index[i]!=-1)
343  {
344  Neighbour_map[std::make_pair(local_el_pt,i)] =
345  finite_element_pt(index[i]);
346  }
347  }
348  }
349  }
350 
351  //Now set up the neighbour information
352  /*const unsigned n_element = this->nelement();
353  for(unsigned e=0;e<n_element;e++)
354  {
355  dynamic_cast<ELEMENT*>(this->element_pt(e))->setup_face_neighbour_info();
356  }*/
357  this->setup_face_neighbour_info();
358  }
Eigen::bfloat16_impl::tan
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 tan(const bfloat16 &a)
Definition: BFloat16.h:633

References e(), i, Global_Parameters::Lx, Global_Parameters::Ly, n, oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), GlobalParameters::Nx, GlobalParameters::Ny, constants::pi, oomph::MathematicalConstants::Pi, Eigen::bfloat16_impl::tan(), BiharmonicTestFunctions2::theta, and oomph::Node::x().

◆ TwoDDGMesh() [5/5]

template<class ELEMENT >
TwoDDGMesh< ELEMENT >::TwoDDGMesh ( const unsigned &  Nx,
const unsigned &  Ny,
TimeStepper *  time_stepper_pt = &Mesh::Default_TimeStepper 
)
inline
91  {
92  Vector<double> s_fraction;
93  //Lengths of the mesh
94  double Lx = 10.0;
95  double Ly = 10.0;
96 
97  double llx = -5.0;
98  double lly = -5.0;
99  //Work out the length of each element in the x and y directions
100  //(Assuming uniform spacing)
101  double el_length[2] = {Lx/(double)Nx, Ly/(double)Ny};
102 
103  //Vector of booleans for boundary information
104  std::vector<bool> boundary_info;
105 
106  //loop over the elements in x
107  for(unsigned ex=0;ex<Nx;ex++)
108  {
109  //loop over the element in y
110  for(unsigned ey=0;ey<Ny;ey++)
111  {
112  //Create a new DG element
113  ELEMENT* local_element_pt = new ELEMENT;
114  //Find the number of nodes
115  const unsigned n_node = local_element_pt->nnode();
116  //Have we constructed
117  bool constructed = false;
118 
119  //If we are on a boundary
120  if(ex==0)
121  {
122  const unsigned n_p = local_element_pt->nnode_1d();
123  boundary_info.resize(n_node);
124  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
125  //The left hand side is on a boundary
126  for(unsigned n=0;n<n_p;n++)
127  {
128  boundary_info[n*n_p] = true;
129  }
130 
131  //Lower left
132  if(ey==0)
133  {
134  for(unsigned n=1;n<n_p;n++) {boundary_info[n] = true;}
135  }
136 
137  //Upper left
138  if(ey==Ny-1)
139  {
140  for(unsigned n=1;n<n_p;n++) {boundary_info[n_p*(n_p-1) + n] = true;}
141  }
142 
143  //Now construct
144  local_element_pt->construct_boundary_nodes_and_faces(
145  this, boundary_info, time_stepper_pt);
146 
147  constructed = true;
148  }
149  //Right boundary
150  else if(ex==Nx-1)
151  {
152  const unsigned n_p = local_element_pt->nnode_1d();
153  boundary_info.resize(n_node);
154  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
155  //The right-hand side is on a boundary
156  for(unsigned n=0;n<n_p;n++)
157  {
158  boundary_info[n_p-1 + n*n_p] = true;
159  }
160 
161  //Lower right
162  if(ey==0)
163  {
164  for(unsigned n=0;n<n_p-1;n++) {boundary_info[n] = true;}
165  }
166 
167  //Upper right
168  if(ey==Ny-1)
169  {
170  for(unsigned n=0;n<n_p-1;n++)
171  {boundary_info[n_p*(n_p-1) + n] = true;}
172  }
173 
174  //Now construct
175  local_element_pt->construct_boundary_nodes_and_faces(
176  this, boundary_info, time_stepper_pt);
177 
178  constructed = true;
179  }
180  //Otherwise it's in the middle
181  else
182  {
183  //If on the bottom or top
184  if((ey==0) || (ey==Ny-1))
185  {
186  const unsigned n_p = local_element_pt->nnode_1d();
187  boundary_info.resize(n_node);
188  for(unsigned n=0;n<n_node;n++) {boundary_info[n] = false;}
189 
190  //If the bottom
191  if(ey==0)
192  {
193  for(unsigned n=0;n<n_p;n++)
194  {
195  boundary_info[n] = true;
196  }
197  }
198 
199  //If on the top
200  if(ey==Ny-1)
201  {
202  for(unsigned n=0;n<n_p;n++)
203  {boundary_info[n_p*(n_p-1) + n] = true;}
204  }
205 
206 
207  //Now construct
208  local_element_pt->construct_boundary_nodes_and_faces(
209  this, boundary_info, time_stepper_pt);
210 
211  constructed = true;
212  }
213  }
214 
215  //If we haven't already constructed
216  if(!constructed)
217  {
218  //Construct nodes and faces
219  local_element_pt->construct_nodes_and_faces(this,time_stepper_pt);
220  }
221 
222  //Find the lower left corner of the element
223  double ll_corner[2] = {llx + ex*el_length[0],lly + ey*el_length[1]};
224  //For each of the nodes set the position
225  for(unsigned n=0;n<n_node;n++)
226  {
227  //Get pointer to the node
228  Node* nod_pt = local_element_pt->node_pt(n);
229  //Get the relative position in local coordinates
230  local_element_pt->local_fraction_of_node(n,s_fraction);
231  //Loop over the coordinates and set the position
232  for(unsigned i=0;i<2;i++)
233  {
234  nod_pt->x(i) = ll_corner[i] + el_length[i]*s_fraction[i];
235  }
236  //Now set the value
237  nod_pt->set_value(0,1.0);
238 
239  //Add each node to the node list
240  Node_pt.push_back(nod_pt);
241  }
242  //Now add the element to the list
243  Element_pt.push_back(local_element_pt);
244  }
245  }
246 
247  //Set up the boundary nodes
248  this->set_nboundary(4);
249 
250  //Now let's add the boundary nodes
251  //Left-hand and right-hand sides
252  for(unsigned e=0;e<Ny;e++)
253  {
254  //Left
255  FiniteElement* elem_pt = this->finite_element_pt(e);
256  unsigned n_p = elem_pt->nnode_1d();
257  for(unsigned n=0;n<n_p;n++)
258  {
259  this->add_boundary_node(3,elem_pt->node_pt(n_p*n));
260  }
261  //Right
262  elem_pt = this->finite_element_pt(Ny*(Nx-1)+e);
263  n_p = elem_pt->nnode_1d();
264  for(unsigned n=0;n<n_p;n++)
265  {
266  this->add_boundary_node(1,elem_pt->node_pt(n_p-1 + n_p*n));
267  }
268  }
269 
270  //Top and bottom
271  for(unsigned e=0;e<Nx;e++)
272  {
273  //Bottom
274  FiniteElement* elem_pt = this->finite_element_pt(Ny*e);
275  unsigned n_p = elem_pt->nnode_1d();
276  for(unsigned n=0;n<n_p;n++)
277  {
278  this->add_boundary_node(0,elem_pt->node_pt(n));
279  }
280 
281  //Top
282  elem_pt = this->finite_element_pt(Ny-1 + Ny*e);
283  n_p = elem_pt->nnode_1d();
284  for(unsigned n=0;n<n_p;n++)
285  {
286  this->add_boundary_node(2,elem_pt->node_pt((n_p-1)*n_p + n));
287  }
288  }
289 
290  //Now loop over all the elements and set up the neighbour map
291  for(unsigned ex=0;ex<Nx;ex++)
292  {
293  for(unsigned ey=0;ey<Ny;ey++)
294  {
295  //Get pointer to the element
296  unsigned element_index = ex*Ny + ey;
297  FiniteElement* local_el_pt = finite_element_pt(element_index);
298 
299  //Storage for indices of neighbours
300  unsigned index[4];
301 
302  //North neighbour (just one up)
303  if(ey < (Ny-1)) {index[0] = element_index + 1;}
304  //If at top, return index at bottom (periodic conditions)
305  else {index[0] = ex*Ny;}
306 
307  //East neighbour (just one across)
308  if(ex < (Nx-1)) {index[1] = element_index + Ny;}
309  //If at side return index at other side (periodic conditions)
310  else {index[1] = ey;}
311 
312  //South neighbour
313  if(ey > 0) {index[2] = element_index - 1;}
314  else {index[2] = element_index + Ny-1;}
315 
316  //West neighbour
317  if(ex > 0) {index[3] = element_index - Ny;}
318  else {index[3] = element_index + (Nx-1)*Ny;}
319 
320  //Now store the details in the mape
321  for(unsigned i=0;i<4;i++)
322  {
323  Neighbour_map[std::make_pair(local_el_pt,i)] =
324  finite_element_pt(index[i]);
325  }
326  }
327  }
328 
329  //Now set up the neighbour information
330  this->setup_face_neighbour_info();
331  }

References e(), i, Global_Parameters::Lx, Global_Parameters::Ly, n, oomph::FiniteElement::nnode_1d(), oomph::FiniteElement::node_pt(), GlobalParameters::Nx, GlobalParameters::Ny, oomph::Data::set_value(), and oomph::Node::x().

Member Function Documentation

◆ neighbour_finder() [1/5]

template<class ELEMENT >
void TwoDDGMesh< ELEMENT >::neighbour_finder ( FiniteElement *const &  bulk_element_pt,
const int &  face_index,
const Vector< double > &  s_bulk,
FaceElement *&  face_element_pt,
Vector< double > &  s_face 
)
inlinevirtual

Reimplemented from oomph::DGMesh.

173  {
174  //We have a single face coordinate of size 1
175  s_face.resize(1);
176 
177  //Now things differ depending upon which face we are located
178  switch(face_index)
179  {
180  //North face
181  case 2:
182  //The neighbouring face element is the south face of the neighbour
183  face_element_pt =
184  dynamic_cast<ELEMENT*>(
185  Neighbour_map[std::make_pair(bulk_element_pt,0)])->
186  face_element_pt(2);
187  //Then set the face coordinate
188  s_face[0] = s_bulk[0];
189  break;
190 
191  //East face
192  case 1:
193  //The neighbouring face element is the west faec of the neighbour
194  face_element_pt =
195  dynamic_cast<ELEMENT*>(
196  Neighbour_map[std::make_pair(bulk_element_pt,1)])->
197  face_element_pt(3);
198  //Then set the face coordinate
199  s_face[0] = s_bulk[1];
200  break;
201 
202  //South face
203  case -2:
204  //the neighbouring face element is the north face of the neighbour
205  face_element_pt =
206  dynamic_cast<ELEMENT*>(
207  Neighbour_map[std::make_pair(bulk_element_pt,2)])
208  ->face_element_pt(0);
209  //Then set the face coordiante
210  s_face[0] = s_bulk[0];
211  break;
212 
213  //West face
214  case -1:
215  //The neighbouring face element is the east face of the neighbour
216  face_element_pt =
217  dynamic_cast<ELEMENT*>
218  (Neighbour_map[std::make_pair(bulk_element_pt,3)])->
219  face_element_pt(1);
220  s_face[0] = s_bulk[1];
221  break;
222 
223  default:
224  throw OomphLibError("Coordinate is on no face or not on a unique face",
225  OOMPH_CURRENT_FUNCTION,
226  OOMPH_EXCEPTION_LOCATION);
227  break;
228  }
229  }
oomph::OomphLibError
Definition: oomph_definitions.h:222
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ neighbour_finder() [2/5]

template<class ELEMENT >
void TwoDDGMesh< ELEMENT >::neighbour_finder ( FiniteElement *const &  bulk_element_pt,
const int &  face_index,
const Vector< double > &  s_bulk,
FaceElement *&  face_element_pt,
Vector< double > &  s_face 
)
inlinevirtual

Reimplemented from oomph::DGMesh.

310  {
311  //We have a single face coordinate of size 1
312  s_face.resize(1);
313 
314  //Now things differ depending upon which face we are located
315  switch(face_index)
316  {
317  //North face
318  case 2:
319  //The neighbouring face element is the south face of the neighbour
320  face_element_pt =
321  dynamic_cast<ELEMENT*>(
322  Neighbour_map[std::make_pair(bulk_element_pt,0)])->
323  face_element_pt(2);
324  //Then set the face coordinate
325  s_face[0] = s_bulk[0];
326  break;
327 
328  //East face
329  case 1:
330  //The neighbouring face element is the west faec of the neighbour
331  face_element_pt =
332  dynamic_cast<ELEMENT*>(
333  Neighbour_map[std::make_pair(bulk_element_pt,1)])->
334  face_element_pt(3);
335  //Then set the face coordinate
336  s_face[0] = s_bulk[1];
337  break;
338 
339  //South face
340  case -2:
341  //the neighbouring face element is the north face of the neighbour
342  face_element_pt =
343  dynamic_cast<ELEMENT*>(
344  Neighbour_map[std::make_pair(bulk_element_pt,2)])
345  ->face_element_pt(0);
346  //Then set the face coordiante
347  s_face[0] = s_bulk[0];
348  break;
349 
350  //West face
351  case -1:
352  //The neighbouring face element is the east face of the neighbour
353  face_element_pt =
354  dynamic_cast<ELEMENT*>
355  (Neighbour_map[std::make_pair(bulk_element_pt,3)])->
356  face_element_pt(1);
357  s_face[0] = s_bulk[1];
358  break;
359 
360  default:
361  throw OomphLibError("Coordinate is on no face or not on a unique face",
362  OOMPH_CURRENT_FUNCTION,
363  OOMPH_EXCEPTION_LOCATION);
364  break;
365  }
366  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ neighbour_finder() [3/5]

template<class ELEMENT >
void TwoDDGMesh< ELEMENT >::neighbour_finder ( FiniteElement *const &  bulk_element_pt,
const int &  face_index,
const Vector< double > &  s_bulk,
FaceElement *&  face_element_pt,
Vector< double > &  s_face 
)
inlinevirtual

Reimplemented from oomph::DGMesh.

687  {
688  //We have a single face coordinate of size 1
689  s_face.resize(1);
690 
691  ELEMENT* elem_pt=0;
692 
693  //Now things differ depending upon which face we are located
694  switch(face_index)
695  {
696  //North face
697  case 2:
698  //Get the neighbour
699  elem_pt = dynamic_cast<ELEMENT*>(
700  Neighbour_map[std::make_pair(bulk_element_pt,0)]);
701 
702  if(elem_pt!=0)
703  {
704  //The neighbouring face element is the south face of the neighbour
705  face_element_pt = elem_pt->face_element_pt(2);
706  }
707  //Otherwise we build in the face reflection element
708  else
709  {
710  /*face_element_pt = new DGEulerFaceReflectionElement<ELEMENT>(
711  bulk_element_pt,2);
712  Face_element_pt.push_back(face_element_pt);*/
713 
714  //dynamic_cast<ELEMENT*>(bulk_element_pt)->face_element_pt(0);
715  face_element_pt = Face_element_pt[std::make_pair(bulk_element_pt,0)];
716  if(face_element_pt==0)
717  {
718  throw OomphLibError("Something wrong",
719  OOMPH_CURRENT_FUNCTION,
720  OOMPH_EXCEPTION_LOCATION);
721  }
722  }
723 
724  //Then set the face coordinate
725  s_face[0] = s_bulk[0];
726  break;
727 
728  //East face
729  case 1:
730  //Get the neighbour
731  elem_pt = dynamic_cast<ELEMENT*>(
732  Neighbour_map[std::make_pair(bulk_element_pt,1)]);
733 
734  //The neighbouring face element is the west faec of the neighbour
735  if(elem_pt!=0)
736  {
737  face_element_pt = elem_pt->face_element_pt(3);
738  }
739  else
740  {
741  ELEMENT* cast_bulk_element_pt = dynamic_cast<ELEMENT*>(bulk_element_pt);
742  //If we're on the outlet set the face to be the face
743  if(cast_bulk_element_pt->face_element_pt(1)->node_pt(0)->
744  is_on_boundary(1))
745  {
746  //std::cout << cast_bulk_element_pt->face_element_pt(1)->node_pt(0)->x(0) << " "
747  // << cast_bulk_element_pt->face_element_pt(1)->node_pt(0)->x(1) << "\n ";
748 
749  face_element_pt =
750  dynamic_cast<ELEMENT*>(bulk_element_pt)->face_element_pt(1);
751  }
752  //Otherwise we must be on the corner
753  else
754  {
755  //face_element_pt = new DGEulerFaceReflectionElement<ELEMENT>(
756  // bulk_element_pt,1);
757  //Face_element_pt.push_back(face_element_pt);
758  face_element_pt = Face_element_pt[std::make_pair(bulk_element_pt,1)];
759  if(face_element_pt==0)
760  {
761  throw OomphLibError("Something wrong",
762  OOMPH_CURRENT_FUNCTION,
763  OOMPH_EXCEPTION_LOCATION);
764  }
765 
766  }
767  }
768 
769  //Then set the face coordinate
770  s_face[0] = s_bulk[1];
771  break;
772 
773  //South face
774  case -2:
775  //Get the neighbour
776  elem_pt = dynamic_cast<ELEMENT*>(
777  Neighbour_map[std::make_pair(bulk_element_pt,2)]);
778 
779  //the neighbouring face element is the north face of the neighbour
780  if(elem_pt!=0)
781  {
782  face_element_pt = elem_pt->face_element_pt(0);
783  }
784  else
785  {
786  //Build in the face reflection element
787  //face_element_pt = new DGEulerFaceReflectionElement<ELEMENT>(
788  // bulk_element_pt,-2);
789  //Face_element_pt.push_back(face_element_pt);
790  //face_element_pt =
791  // dynamic_cast<ELEMENT*>(bulk_element_pt)->face_element_pt(2);
792 
793 
794  face_element_pt = Face_element_pt[std::make_pair(bulk_element_pt,2)];
795 
796  if(face_element_pt==0)
797  {
798  throw OomphLibError("Something wrong",
799  OOMPH_CURRENT_FUNCTION,
800  OOMPH_EXCEPTION_LOCATION);
801  }
802 
803  }
804 
805  //Then set the face coordiante
806  s_face[0] = s_bulk[0];
807  break;
808 
809  //West face
810  case -1:
811  //Get the neighbour
812  elem_pt = dynamic_cast<ELEMENT*>(
813  Neighbour_map[std::make_pair(bulk_element_pt,3)]);
814 
815  //The neighbouring face element is the east face of the neighbour
816  if(elem_pt!=0)
817  {
818  face_element_pt = elem_pt->face_element_pt(1);
819  }
820  else
821  {
822  face_element_pt =
823  dynamic_cast<ELEMENT*>(bulk_element_pt)->face_element_pt(3);
824  }
825  s_face[0] = s_bulk[1];
826  break;
827 
828  default:
829  throw OomphLibError("Coordinate is on no face or not on a unique face",
830  OOMPH_CURRENT_FUNCTION,
831  OOMPH_EXCEPTION_LOCATION);
832  break;
833  }
834  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ neighbour_finder() [4/5]

template<class ELEMENT >
void TwoDDGMesh< ELEMENT >::neighbour_finder ( FiniteElement *const &  bulk_element_pt,
const int &  face_index,
const Vector< double > &  s_bulk,
FaceElement *&  face_element_pt,
Vector< double > &  s_face 
)
inlinevirtual

Reimplemented from oomph::DGMesh.

366  {
367  //We have a single face coordinate of size 1
368  s_face.resize(1);
369 
370  ELEMENT* elem_pt = 0;
371 
372  //Now things differ depending upon which face we are located
373  switch(face_index)
374  {
375  //North face
376  case 2:
377  //Get the neighbour
378  elem_pt = dynamic_cast<ELEMENT*>(
379  Neighbour_map[std::make_pair(bulk_element_pt,0)]);
380 
381  if(elem_pt!=0)
382  {
383  //The neighbouring face element is the south face of the neighbour
384  face_element_pt = elem_pt->face_element_pt(2);
385  }
386  //Otherwise free outflow
387  else
388  {
389  face_element_pt =
390  dynamic_cast<ELEMENT*>(bulk_element_pt)->face_element_pt(0);
391  }
392 
393  //Then set the face coordinate
394  s_face[0] = s_bulk[0];
395  break;
396 
397  //East face
398  case 1:
399  //Get the neighbour
400  elem_pt =
401  dynamic_cast<ELEMENT*>(
402  Neighbour_map[std::make_pair(bulk_element_pt,1)]);
403 
404  if(elem_pt!=0)
405  {
406  face_element_pt = elem_pt->face_element_pt(3);
407  }
408  else
409  {
410  face_element_pt =
411  dynamic_cast<ELEMENT*>(bulk_element_pt)->face_element_pt(1);
412  }
413 
414  //Then set the face coordinate
415  s_face[0] = s_bulk[1];
416  break;
417 
418  //South face
419  case -2:
420  //Get the neighbour
421  elem_pt =
422  dynamic_cast<ELEMENT*>(
423  Neighbour_map[std::make_pair(bulk_element_pt,2)]);
424 
425  if(elem_pt!=0)
426  {
427  face_element_pt = elem_pt->face_element_pt(0);
428  }
429  //Otherwise use reflection element
430  else
431  {
432  face_element_pt = Face_element_pt[std::make_pair(bulk_element_pt,2)];
433  if(face_element_pt==0)
434  {
435  throw OomphLibError("Something wrong",
436  OOMPH_CURRENT_FUNCTION,
437  OOMPH_EXCEPTION_LOCATION);
438  }
439  }
440 
441  //Then set the face coordiante
442  s_face[0] = s_bulk[0];
443  break;
444 
445  //West face
446  case -1:
447  //Get the neighbour
448  elem_pt =
449  dynamic_cast<ELEMENT*>
450  (Neighbour_map[std::make_pair(bulk_element_pt,3)]);
451 
452  if(elem_pt!=0)
453  {
454  face_element_pt = elem_pt->face_element_pt(1);
455  }
456  else
457  {
458  face_element_pt = dynamic_cast<ELEMENT*>(bulk_element_pt)
459  ->face_element_pt(3);
460  }
461 
462  s_face[0] = s_bulk[1];
463  break;
464 
465  default:
466  throw OomphLibError("Coordinate is on no face or not on a unique face",
467  OOMPH_CURRENT_FUNCTION,
468  OOMPH_EXCEPTION_LOCATION);
469  break;
470  }
471  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

◆ neighbour_finder() [5/5]

template<class ELEMENT >
void TwoDDGMesh< ELEMENT >::neighbour_finder ( FiniteElement *const &  bulk_element_pt,
const int &  face_index,
const Vector< double > &  s_bulk,
FaceElement *&  face_element_pt,
Vector< double > &  s_face 
)
inlinevirtual

Reimplemented from oomph::DGMesh.

339  {
340  //We have a single face coordinate of size 1
341  s_face.resize(1);
342 
343  //Now things differ depending upon which face we are located
344  switch(face_index)
345  {
346  //North face
347  case 2:
348  //The neighbouring face element is the south face of the neighbour
349  face_element_pt =
350  dynamic_cast<ELEMENT*>(
351  Neighbour_map[std::make_pair(bulk_element_pt,0)])->
352  face_element_pt(2);
353  //Then set the face coordinate
354  s_face[0] = s_bulk[0];
355  break;
356 
357  //East face
358  case 1:
359  //The neighbouring face element is the west faec of the neighbour
360  face_element_pt =
361  dynamic_cast<ELEMENT*>(
362  Neighbour_map[std::make_pair(bulk_element_pt,1)])->
363  face_element_pt(3);
364  //Then set the face coordinate
365  s_face[0] = s_bulk[1];
366  break;
367 
368  //South face
369  case -2:
370  //the neighbouring face element is the north face of the neighbour
371  face_element_pt =
372  dynamic_cast<ELEMENT*>(
373  Neighbour_map[std::make_pair(bulk_element_pt,2)])
374  ->face_element_pt(0);
375  //Then set the face coordiante
376  s_face[0] = s_bulk[0];
377  break;
378 
379  //West face
380  case -1:
381  //The neighbouring face element is the east face of the neighbour
382  face_element_pt =
383  dynamic_cast<ELEMENT*>
384  (Neighbour_map[std::make_pair(bulk_element_pt,3)])->
385  face_element_pt(1);
386  s_face[0] = s_bulk[1];
387  break;
388 
389  default:
390  throw OomphLibError("Coordinate is on no face or not on a unique face",
391  OOMPH_CURRENT_FUNCTION,
392  OOMPH_EXCEPTION_LOCATION);
393  break;
394  }
395  }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

Member Data Documentation

◆ Face_element_pt

template<class ELEMENT >
std::map< std::pair< FiniteElement *, int >, FaceElement * > TwoDDGMesh< ELEMENT >::Face_element_pt
private

◆ Neighbour_map [1/2]

template<class ELEMENT >
std::map< std::pair< FiniteElement *, unsigned >, FiniteElement * > TwoDDGMesh< ELEMENT >::Neighbour_map
private

◆ Neighbour_map [2/2]

template<class ELEMENT >
std::map<std::pair<FiniteElement*,int>,FiniteElement*> TwoDDGMesh< ELEMENT >::Neighbour_map
private
The documentation for this class was generated from the following files: