CppNoddy  0.92
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | Protected Types | Protected Member Functions | Protected Attributes | List of all members
CppNoddy::TwoD_TVDLF_Mesh Class Reference

#include <TwoD_TVDLF_Mesh.h>

Public Member Functions

 TwoD_TVDLF_Mesh (const DenseVector< double > &X, const DenseVector< double > &Y, TwoD_Hyperbolic_System *ptr, fn_ptr init_ptr)
 Constructor for the Finite Volume Mesh using linear elements. More...
 
virtual ~TwoD_TVDLF_Mesh ()
 Empty desctructor. More...
 
void dump_gnu (std::string filename)
 Dump the data to a given filename in a gnuplot format. More...
 
void dump_nodes_x (std::string filename) const
 Dump the x-nodal positions to a given filename. More...
 
void dump_nodes_y (std::string filename) const
 Dump the y-nodal positions to a given filename. More...
 
void dump_data (std::string filename)
 Dump the data over the nodal positions to a given filename. More...
 
void set_limiter (const unsigned &id)
 Set the limiter type to be applied in the slope values. More...
 
double update (const double &CFL, const double &max_dt=std::numeric_limits< long double >::max())
 Update the mesh object. More...
 
void update_to (const double &CFL, const double &t_end)
 Update the mesh object to a set time level. More...
 
DenseVector< double > integrate (std::string mesh_colour="black")
 Integrate the concentration values across the entire mesh. More...
 
elt_iter get_elt_iter_from_x (const DenseVector< double > &x, std::string mesh_colour="black")
 Given a global coordinate, return a pointer to the elt that contains that point. More...
 
DenseVector< double > get_point_values (const DenseVector< double > &x)
 Get the vector of unknowns at a point in the 2D mesh. More...
 
const double & get_time () const
 Get a const reference to the time value for the current mesh. More...
 
virtual void actions_before_time_step1 (const double &time_step)
 A virtual method that is run before the first time update. More...
 
virtual void actions_before_time_step2 (const double &time_step)
 A virtual method that is run before the second time update. More...
 

Public Attributes

vector_of_elts BLACK_ELTS
 An STL vector of linear elements – the black mesh. More...
 
vector_of_elts RED_ELTS
 An STL vector of linear elements – the red mesh. More...
 

Protected Types

typedef std::vector< TwoD_TVDLF_Elt > vector_of_elts
 iterators for the vector of elements More...
 
typedef vector_of_elts::const_iterator celt_iter
 
typedef vector_of_elts::iterator elt_iter
 
typedef std::vector< TwoD_TVDLF_Elt * > vector_of_boundary_elts
 
typedef vector_of_boundary_elts::iterator bdry_elt_iter
 
typedef void(* fn_ptr) (const double &, const double &, DenseVector< double > &)
 function pointer used in the initial conditions More...
 

Protected Member Functions

elt_iter get_elt_iter_from_elt_iter (elt_iter e, std::string target_colour, int corner_index)
 Given an element in the INITIAL STRUCTURED MESH, this method will return an iterator to an element in the other mesh that overlaps the corner specified by corner_index. More...
 
vector_of_elts * get_elts_from_colour (std::string mesh_colour)
 Given a choice of black or red mesh, return a pointer to the appropriate vector of elements. More...
 
std::size_t get_number_elts_in_x (std::string mesh_colour)
 Given a choice of black or red mesh, return the number of elements in the x-direction. More...
 
void calc_slopes (vector_of_elts *elt_vector)
 Use the appropriate limiter to approximate the slope in each element in the mesh. More...
 
DenseVector< double > east_diff (elt_iter e) const
 Compute a finite difference approximation of the derivative in the compass direction. More...
 
DenseVector< double > west_diff (elt_iter e) const
 Compute a finite difference approximation of the derivative in the compass direction. More...
 
DenseVector< double > north_diff (elt_iter e) const
 Compute a finite difference approximation of the derivative in the compass direction. More...
 
DenseVector< double > south_diff (elt_iter e) const
 Compute a finite difference approximation of the derivative in the compass direction. More...
 
DenseVector< double > NS_diff (elt_iter e)
 Compute a finite difference approximation of the derivative in the compass direction. More...
 
DenseVector< double > EW_diff (elt_iter e)
 Compute a finite difference approximation of the derivative in the compass direction. More...
 
int sgn (double a) const
 Sign of a double. More...
 
DenseVector< double > minmod (DenseVector< double > A, DenseVector< double > B) const
 A vector version of the minmod operator. More...
 
DenseVector< double > maxmod (DenseVector< double > A, DenseVector< double > B) const
 A vector version of the maxmod operator. More...
 
void boundary_diff (elt_iter e, const int &face_index, DenseVector< double > &diff) const
 Given an element iterator and the local coordinate this will return zero for any components specified as inflow boundary conditions in line with Levy & Tadmor (1997) and set the centre nodal value to the edge value. More...
 
void set_boundary_Q (vector_of_elts *elts)
 Loops over all boundary elements and sets the Q values in each one to be the value specified by the edge_values method. More...
 

Protected Attributes

unsigned LIMITER
 Slope limiter method. More...
 
std::size_t ORDER_OF_SYSTEM
 order of the conservative system More...
 
std::size_t NX
 number of faces in the x & y directions More...
 
std::size_t NY
 
double MESH_TIME
 the time level of the mesh More...
 
fn_ptr p_Q_INIT
 function pointer to a funnction that defines the initial distribution More...
 

Detailed Description

Definition at line 19 of file TwoD_TVDLF_Mesh.h.

Member Typedef Documentation

◆ bdry_elt_iter

typedef vector_of_boundary_elts::iterator CppNoddy::TwoD_TVDLF_Mesh::bdry_elt_iter
protected

Definition at line 27 of file TwoD_TVDLF_Mesh.h.

◆ celt_iter

typedef vector_of_elts::const_iterator CppNoddy::TwoD_TVDLF_Mesh::celt_iter
protected

Definition at line 23 of file TwoD_TVDLF_Mesh.h.

◆ elt_iter

typedef vector_of_elts::iterator CppNoddy::TwoD_TVDLF_Mesh::elt_iter
protected

Definition at line 24 of file TwoD_TVDLF_Mesh.h.

◆ fn_ptr

typedef void(* CppNoddy::TwoD_TVDLF_Mesh::fn_ptr) (const double &, const double &, DenseVector< double > &)
protected

function pointer used in the initial conditions

Definition at line 30 of file TwoD_TVDLF_Mesh.h.

◆ vector_of_boundary_elts

typedef std::vector<TwoD_TVDLF_Elt*> CppNoddy::TwoD_TVDLF_Mesh::vector_of_boundary_elts
protected

Definition at line 26 of file TwoD_TVDLF_Mesh.h.

◆ vector_of_elts

typedef std::vector<TwoD_TVDLF_Elt> CppNoddy::TwoD_TVDLF_Mesh::vector_of_elts
protected

iterators for the vector of elements

Definition at line 22 of file TwoD_TVDLF_Mesh.h.

Constructor & Destructor Documentation

◆ TwoD_TVDLF_Mesh()

CppNoddy::TwoD_TVDLF_Mesh::TwoD_TVDLF_Mesh ( const DenseVector< double > &  X,
const DenseVector< double > &  Y,
TwoD_Hyperbolic_System *  ptr,
fn_ptr  init_ptr 
)

Constructor for the Finite Volume Mesh using linear elements.

Parameters
XA vector of nodal locations at which the element FACES will positioned
YA vector of nodal locations at which the element FACES will positioned
ptrA pointer to the hyperbolic system applied to this mesh
init_ptrA pointer to a function that defines the initial conditions

Definition at line 16 of file TwoD_TVDLF_Mesh.cpp.

18 {
19
20#ifdef DEBUG
21 std::cout << "DEBUG: Starting construction of a TwoD_TVDLF_Mesh object. \n";
22#endif
23 MESH_TIME = 0.0;
24 NX = X.size();
25 NY = Y.size();
26 if(std::min(NX - 1, NY - 1) <= 1) {
27 std::string problem;
28 problem = " The TwoD_TVDLF_Mesh object is trying to construct itself \n";
29 problem += " with just one element in one of the directions! \n";
30 throw ExceptionRuntime(problem);
31 }
32
33#ifdef DEBUG
34 std::cout << "DEBUG: configuration of the black mesh \n";
35#endif
36
37 // reserve appropriate space to avoid re-allocation later
38 BLACK_ELTS.reserve((NX - 1) * (NY - 1));
39 RED_ELTS.reserve(NX * NY);
40
41 // set up the fn ptr to the initial conditions fn
42 p_Q_INIT = init_ptr;
43 // default limiter
44 LIMITER = 0;
45 // store the order of the conservative system here for simplicity
46 ORDER_OF_SYSTEM = ptr -> get_order();
47
48 // face index of edges
49 std::set<int> faces_s;
50 faces_s.insert(0);
51 std::set<int> faces_e;
52 faces_e.insert(1);
53 std::set<int> faces_n;
54 faces_n.insert(2);
55 std::set<int> faces_w;
56 faces_w.insert(3);
57 // face indices of corners
58 std::set<int> faces_se(faces_s);
59 faces_se.insert(1);
60 std::set<int> faces_ne(faces_n);
61 faces_ne.insert(1);
62 std::set<int> faces_nw(faces_n);
63 faces_nw.insert(3);
64 std::set<int> faces_sw(faces_s);
65 faces_sw.insert(3);
66 // set up the black elements
67 {
68 // southern row
69 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[0], X[1], Y[0], Y[1], ptr, true, faces_sw));
70 for(std::size_t i = 1; i <= NX - 3; ++i) {
71 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[i], X[i+1], Y[0], Y[1], ptr, true, faces_s));
72 }
73 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[NX-2], X[NX-1], Y[0], Y[1], ptr, true, faces_se));
74 // interior elts
75 for(std::size_t j = 1; j <= NY - 3; ++j) {
76 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[0], X[1], Y[j], Y[j+1], ptr, true, faces_w));
77 for(std::size_t i = 1; i <= NX - 3; ++i) {
78 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[i], X[i+1], Y[j], Y[j+1], ptr));
79 }
80 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[NX-2], X[NX-1], Y[j], Y[j+1], ptr, true, faces_e));
81 }
82 // northern row
83 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[0], X[1], Y[NY-2], Y[NY-1], ptr, true, faces_nw));
84 for(std::size_t i = 1; i <= NX - 3; ++i) {
85 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[i], X[i+1], Y[NY-2], Y[NY-1], ptr, true, faces_n));
86 }
87 BLACK_ELTS.push_back(TwoD_TVDLF_Elt(X[NX-2], X[NX-1], Y[NY-2], Y[NY-1], ptr, true, faces_ne));
88 }
89
90#ifdef DEBUG
91 std::cout << "DEBUG: configuration of the red mesh \n";
92#endif
93
94 // set up the red elements
95 {
96 // southern row
97 RED_ELTS.push_back(TwoD_TVDLF_Elt(X[0], (X[1] + X[0]) / 2, Y[0], (Y[0] + Y[1]) / 2, ptr, true, faces_sw));
98 for(std::size_t i = 1; i <= NX - 2; ++i) {
99 RED_ELTS.push_back(TwoD_TVDLF_Elt((X[i-1] + X[i]) / 2, (X[i] + X[i+1]) / 2, Y[0], (Y[0] + Y[1]) / 2, ptr, true, faces_s));
100 }
101 RED_ELTS.push_back(TwoD_TVDLF_Elt((X[NX-2] + X[NX-1]) / 2, X[NX-1], Y[0], (Y[0] + Y[1]) / 2, ptr, true, faces_se));
102 // interior elts
103 for(std::size_t j = 1; j <= NY - 2; ++j) {
104 RED_ELTS.push_back(TwoD_TVDLF_Elt(X[0], (X[1] + X[0]) / 2, (Y[j-1] + Y[j]) / 2, (Y[j] + Y[j+1]) / 2, ptr, true, faces_w));
105 for(std::size_t i = 1; i <= NX - 2; ++i) {
106 RED_ELTS.push_back(TwoD_TVDLF_Elt((X[i-1] + X[i]) / 2, (X[i] + X[i+1]) / 2, (Y[j-1] + Y[j]) / 2, (Y[j] + Y[j+1]) / 2, ptr));
107 }
108 RED_ELTS.push_back(TwoD_TVDLF_Elt((X[NX-2] + X[NX-1]) / 2, X[NX-1], (Y[j-1] + Y[j]) / 2, (Y[j] + Y[j+1]) / 2, ptr, true, faces_e));
109 }
110 // northern row
111 RED_ELTS.push_back(TwoD_TVDLF_Elt(X[0], (X[1] + X[0]) / 2, (Y[NY-2] + Y[NY-1]) / 2, Y[NY-1], ptr, true, faces_nw));
112 for(std::size_t i = 1; i <= NX - 2; ++i) {
113 RED_ELTS.push_back(TwoD_TVDLF_Elt((X[i-1] + X[i]) / 2, (X[i] + X[i+1]) / 2, (Y[NY-2] + Y[NY-1]) / 2, Y[NY-1], ptr, true, faces_n));
114 }
115 RED_ELTS.push_back(TwoD_TVDLF_Elt((X[NX-2] + X[NX-1]) / 2, X[NX-1], (Y[NY-2] + Y[NY-1]) / 2, Y[NY-1], ptr, true, faces_ne));
116 }
117
118#ifdef DEBUG
119 std::cout << "DEBUG: computing black to red projection \n";
120#endif
121
122 // Now we need to pre-compute the mesh projection from black-red-black.
123 // Each element needs to know which element in the other mesh contributes
124 // to it in the projection scheme. It also needs to know which of its faces
125 // the element contributes to in the flux computation.
126 // We assume the constructor is not required to be efficient & just brute
127 // force it here rather than working out all the index mappings by hand.
128 //
129 // loop through all the black elts
130 elt_iter eb = BLACK_ELTS.begin();
131 while(eb != BLACK_ELTS.end()) {
132 // local coordinates of the corners
133 DenseVector<double> ne(2, 1.0);
134 DenseVector<double> sw(-ne);
135 DenseVector<double> se(2, 1.0);
136 se[ 1 ] = -1.0;
137 DenseVector<double> nw(-se);
138 // global coordinates of the corners
139 DenseVector<double> bx_sw = eb -> get_x(sw);
140 DenseVector<double> bx_ne = eb -> get_x(ne);
141 DenseVector<double> bx_nw = eb -> get_x(nw);
142 DenseVector<double> bx_se = eb -> get_x(se);
143 // get an iterator to the red elt that contains each corner
144 // & add it as a contribution
145 //
146 // nw corner
147 {
148 //elt_iter er = get_elt_iter_from_x( bx_nw, "red" );
149 elt_iter er = get_elt_iter_from_elt_iter(eb, "red", 3);
150 DenseVector<double> s = er -> get_s(bx_nw);
151 DenseVector<double> s1(2, 0.0);
152 s1[ 0 ] = s[ 0 ];
153 s1[ 1 ] = -1.0;
154 DenseVector<double> s2(2, 0.0);
155 s2[ 0 ] = 1.0;
156 s2[ 1 ] = s[ 1 ];
157 eb -> add_contribution(&(*er), s1, s2, faces_nw);
158 }
159 // sw corner
160 {
161 //elt_iter er = get_elt_iter_from_x( bx_sw, "red" );
162 elt_iter er = get_elt_iter_from_elt_iter(eb, "red", 0);
163 DenseVector<double> s = er -> get_s(bx_sw);
164 DenseVector<double> s1(2, 0.0);
165 s1[ 0 ] = s[ 0 ];
166 s1[ 1 ] = s[ 1 ];
167 DenseVector<double> s2(2, 0.0);
168 s2[ 0 ] = 1.0;
169 s2[ 1 ] = 1.0;
170 eb -> add_contribution(&(*er), s1, s2, faces_sw);
171 }
172 // ne corner
173 {
174 //elt_iter er = get_elt_iter_from_x( bx_ne, "red" );
175 elt_iter er = get_elt_iter_from_elt_iter(eb, "red", 2);
176 DenseVector<double> s = er -> get_s(bx_ne);
177 DenseVector<double> s1(2, 0.0);
178 s1[ 0 ] = -1.0;
179 s1[ 1 ] = -1.0;
180 DenseVector<double> s2(2, 0.0);
181 s2[ 0 ] = s[ 0 ];
182 s2[ 1 ] = s[ 1 ];
183 eb -> add_contribution(&(*er), s1, s2, faces_ne);
184 }
185 // se corner
186 {
187 //elt_iter er = get_elt_iter_from_x( bx_se, "red" );
188 elt_iter er = get_elt_iter_from_elt_iter(eb, "red", 1);
189 DenseVector<double> s = er -> get_s(bx_se);
190 DenseVector<double> s1(2, 0.0);
191 s1[ 0 ] = -1.0;
192 s1[ 1 ] = s[ 1 ];
193 DenseVector<double> s2(2, 0.0);
194 s2[ 0 ] = s[ 0 ];
195 s2[ 1 ] = 1.0;
196 eb -> add_contribution(&(*er), s1, s2, faces_se);
197 }
198 ++eb;
199 }
200
201#ifdef DEBUG
202 std::cout << "DEBUG: computing red to black projection \n";
203#endif
204
205 // now we need to set up the projection from red to black elts.
206 // it's a little more tricky here as we have to avoid the external nodes.
207 elt_iter er = RED_ELTS.begin();
208 while(er != RED_ELTS.end()) {
209 // local coordinates of the corners
210 DenseVector<double> ne(2, 1.0);
211 DenseVector<double> sw(-ne);
212 DenseVector<double> se(2, 1.0);
213 se[ 1 ] = -1.0;
214 DenseVector<double> nw(-se);
215 // global coordinates of the corners
216 DenseVector<double> bx_sw = er -> get_x(sw);
217 DenseVector<double> bx_ne = er -> get_x(ne);
218 DenseVector<double> bx_nw = er -> get_x(nw);
219 DenseVector<double> bx_se = er -> get_x(se);
220 //
221 std::set<int> external = er -> get_external_faces();
222 if((external.find(2) == external.end()) &&
223 (external.find(3) == external.end())) {
224 // nw corner is internal
225 //elt_iter eb = get_elt_iter_from_x( bx_nw, "black" );
226 elt_iter eb = get_elt_iter_from_elt_iter(er, "black", 3);
227 DenseVector<double> s = eb -> get_s(bx_nw);
228 DenseVector<double> s1(2, 0.0);
229 s1[ 0 ] = s[ 0 ];
230 s1[ 1 ] = -1.0;
231 DenseVector<double> s2(2, 0.0);
232 s2[ 0 ] = 1.0;
233 s2[ 1 ] = s[ 1 ];
234 er -> add_contribution(&(*eb), s1, s2, faces_nw);
235 }
236 if((external.find(0) == external.end()) &&
237 (external.find(3) == external.end())) {
238 // sw corner is internal
239 //elt_iter eb = get_elt_iter_from_x( bx_sw, "black" );
240 elt_iter eb = get_elt_iter_from_elt_iter(er, "black", 0);
241 DenseVector<double> s = eb -> get_s(bx_sw);
242 DenseVector<double> s1(2, 0.0);
243 s1[ 0 ] = s[ 0 ];
244 s1[ 1 ] = s[ 1 ];
245 DenseVector<double> s2(2, 0.0);
246 s2[ 0 ] = 1.0;
247 s2[ 1 ] = 1.0;
248 er -> add_contribution(&(*eb), s1, s2, faces_sw);
249 }
250 if((external.find(1) == external.end()) &&
251 (external.find(2) == external.end())) {
252 // ne corner is internal
253 //elt_iter eb = get_elt_iter_from_x( bx_ne, "black" );
254 elt_iter eb = get_elt_iter_from_elt_iter(er, "black", 2);
255 DenseVector<double> s = eb -> get_s(bx_ne);
256 DenseVector<double> s1(2, 0.0);
257 s1[ 0 ] = -1.0;
258 s1[ 1 ] = -1.0;
259 DenseVector<double> s2(2, 0.0);
260 s2[ 0 ] = s[ 0 ];
261 s2[ 1 ] = s[ 1 ];
262 er -> add_contribution(&(*eb), s1, s2, faces_ne);
263 }
264 if((external.find(0) == external.end()) &&
265 (external.find(1) == external.end())) {
266 // se corner is internal
267 //elt_iter eb = get_elt_iter_from_x( bx_se, "black" );
268 elt_iter eb = get_elt_iter_from_elt_iter(er, "black", 1);
269 DenseVector<double> s = eb -> get_s(bx_se);
270 DenseVector<double> s1(2, 0.0);
271 s1[ 0 ] = -1.0;
272 s1[ 1 ] = s[ 1 ];
273 DenseVector<double> s2(2, 0.0);
274 s2[ 0 ] = s[ 0 ];
275 s2[ 1 ] = 1.0;
276 er -> add_contribution(&(*eb), s1, s2, faces_se);
277 }
278 ++er;
279 }
280
281#ifdef DEBUG
282 std::cout << "DEBUG: setting up pointers to NESW elts in black mesh\n";
283#endif
284
285 // we need to set the pointers to neighbouring elts in both meshes
286 eb = BLACK_ELTS.begin();
287 while(eb != BLACK_ELTS.end()) {
288 // the offset for computing N & S elts is no. of faces minus one
289 const std::size_t offset(NX - 1);
290 // get a set of external faces
291 std::set< int > faces(eb -> get_external_faces());
292 if(eb -> face_is_internal(0)) {
293 // southern face is internal
294 eb -> set_ptrs(0, &(*eb) - offset);
295 }
296 if(eb -> face_is_internal(1)) {
297 // eastern face is internal
298 eb -> set_ptrs(1, &(*eb) + 1);
299 }
300 if(eb -> face_is_internal(2)) {
301 // northern face is internal
302 eb -> set_ptrs(2, &(*eb) + offset);
303 }
304 if(eb -> face_is_internal(3)) {
305 // western face is internal
306 eb -> set_ptrs(3, &(*eb) - 1);
307 }
308 ++eb;
309 }
310
311#ifdef DEBUG
312 std::cout << "DEBUG: setting up pointers to NESW elts in red mesh\n";
313#endif
314
315 // we need to set the pointers to neighbouring elts in both meshes
316 er = RED_ELTS.begin();
317 while(er != RED_ELTS.end()) {
318 // the offset for computing N & S elts is no. of faces
319 // as the red mesh has an extra elt per row
320 const std::size_t offset(NX);
321 // get a set of external faces
322 std::set< int > faces(er -> get_external_faces());
323 if(er -> face_is_internal(0)) {
324 // southern face is internal
325 er -> set_ptrs(0, &(*er) - offset);
326 }
327 if(er -> face_is_internal(1)) {
328 // eastern face is internal
329 er -> set_ptrs(1, &(*er) + 1);
330 }
331 if(er -> face_is_internal(2)) {
332 // northern face is internal
333 er -> set_ptrs(2, &(*er) + offset);
334 }
335 if(er -> face_is_internal(3)) {
336 // western face is internal
337 er -> set_ptrs(3, &(*er) - 1);
338 }
339 ++er;
340 }
341
342#ifdef DEBUG
343 std::cout << "DEBUG: initialising the black mesh \n";
344#endif
345
346 // now all that remains is to initialise the starting (black) mesh
347 eb = BLACK_ELTS.begin();
348 while(eb != BLACK_ELTS.end()) {
349 DenseVector<double> s(2, 0.0);
350 // compute to the east
351 s[ 0 ] = 1.0;
352 s[ 1 ] = 0.0;
353 DenseVector<double> xe(eb -> get_x(s));
354 DenseVector<double> Qe(ORDER_OF_SYSTEM, 0.0);
355 p_Q_INIT(xe[0], xe[1], Qe);
356 // compute to the west
357 s[ 0 ] = -1.0;
358 s[ 1 ] = 0.0;
359 DenseVector<double> xw(eb -> get_x(s));
360 DenseVector<double> Qw(ORDER_OF_SYSTEM, 0.0);
361 p_Q_INIT(xw[0], xw[1], Qw);
362 // compute to the north
363 s[ 0 ] = 0.0;
364 s[ 1 ] = 1.0;
365 DenseVector<double> xn(eb -> get_x(s));
366 DenseVector<double> Qn(ORDER_OF_SYSTEM, 0.0);
367 p_Q_INIT(xn[0], xn[1], Qn);
368 // compute to the north
369 s[ 0 ] = 0.0;
370 s[ 1 ] = -1.0;
371 DenseVector<double> xs(eb -> get_x(s));
372 DenseVector<double> Qs(ORDER_OF_SYSTEM, 0.0);
373 p_Q_INIT(xs[0], xs[1], Qs);
374 // difference for the slopes
375 eb -> set_slope_x((Qe - Qw) / (xe[0] - xw[0]));
376 eb -> set_slope_y((Qn - Qs) / (xn[1] - xs[1]));
377 // set the mid value
378 eb -> set_Q_mid((Qe + Qw + Qn + Qs) / 4);
379 ++eb;
380 }
381
382#ifdef DEBUG
383 std::cout << "DEBUG: mesh constructor complete \n";
384#endif
385
386 // now all that remains is to initialise the starting (black) mesh
387 eb = BLACK_ELTS.begin();
388 while(eb != BLACK_ELTS.end()) {
389 DenseVector<double> x(2, 0.0);
390 DenseVector<double> s(2, 0.0);
391 x = eb -> get_x(s);
392 DenseVector<double> Q(ORDER_OF_SYSTEM, 0.0);
393 p_Q_INIT(x[0], x[1], Q);
394 eb -> set_Q_mid(Q);
395 ++eb;
396 }
397
398 calc_slopes(&BLACK_ELTS);
399
400 }
CppNoddy::DenseVector::size
std::size_t size() const
A pass-thru definition to get the size of the vector.
Definition: DenseVector.h:330
CppNoddy::TwoD_TVDLF_Mesh::RED_ELTS
vector_of_elts RED_ELTS
An STL vector of linear elements – the red mesh.
Definition: TwoD_TVDLF_Mesh.h:128
CppNoddy::TwoD_TVDLF_Mesh::calc_slopes
void calc_slopes(vector_of_elts *elt_vector)
Use the appropriate limiter to approximate the slope in each element in the mesh.
Definition: TwoD_TVDLF_Mesh.cpp:690
CppNoddy::TwoD_TVDLF_Mesh::LIMITER
unsigned LIMITER
Slope limiter method.
Definition: TwoD_TVDLF_Mesh.h:373
CppNoddy::TwoD_TVDLF_Mesh::BLACK_ELTS
vector_of_elts BLACK_ELTS
An STL vector of linear elements – the black mesh.
Definition: TwoD_TVDLF_Mesh.h:126
CppNoddy::TwoD_TVDLF_Mesh::get_elt_iter_from_elt_iter
elt_iter get_elt_iter_from_elt_iter(elt_iter e, std::string target_colour, int corner_index)
Given an element in the INITIAL STRUCTURED MESH, this method will return an iterator to an element in...
Definition: TwoD_TVDLF_Mesh.h:139
CppNoddy::TwoD_TVDLF_Mesh::p_Q_INIT
fn_ptr p_Q_INIT
function pointer to a funnction that defines the initial distribution
Definition: TwoD_TVDLF_Mesh.h:381
CppNoddy::TwoD_TVDLF_Mesh::ORDER_OF_SYSTEM
std::size_t ORDER_OF_SYSTEM
order of the conservative system
Definition: TwoD_TVDLF_Mesh.h:375
CppNoddy::TwoD_TVDLF_Mesh::NY
std::size_t NY
Definition: TwoD_TVDLF_Mesh.h:377
CppNoddy::TwoD_TVDLF_Mesh::NX
std::size_t NX
number of faces in the x & y directions
Definition: TwoD_TVDLF_Mesh.h:377
CppNoddy::TwoD_TVDLF_Mesh::MESH_TIME
double MESH_TIME
the time level of the mesh
Definition: TwoD_TVDLF_Mesh.h:379
CppNoddy::TwoD_TVDLF_Mesh::elt_iter
vector_of_elts::iterator elt_iter
Definition: TwoD_TVDLF_Mesh.h:24
CppNoddy::Example::s
double s
relative rotation rate
Definition: BVPKarmanArclength.cpp:30

References BLACK_ELTS, calc_slopes(), get_elt_iter_from_elt_iter(), LIMITER, MESH_TIME, NX, NY, ORDER_OF_SYSTEM, p_Q_INIT, RED_ELTS, and CppNoddy::DenseVector< _Type >::size().

◆ ~TwoD_TVDLF_Mesh()

CppNoddy::TwoD_TVDLF_Mesh::~TwoD_TVDLF_Mesh ( )
virtual

Empty desctructor.

Definition at line 402 of file TwoD_TVDLF_Mesh.cpp.

403 {}

Member Function Documentation

◆ actions_before_time_step1()

virtual void CppNoddy::TwoD_TVDLF_Mesh::actions_before_time_step1 ( const double &  time_step)
inlinevirtual

A virtual method that is run before the first time update.

For user-custom problems.

Parameters
time_stepThe time step size that is about to be taken

Definition at line 110 of file TwoD_TVDLF_Mesh.h.

110 {
111 // Empty by default. If you want to take actions before the time
112 // step then you need to inherit from this basic mesh and implement
113 // this method.
114 }

Referenced by update().

◆ actions_before_time_step2()

virtual void CppNoddy::TwoD_TVDLF_Mesh::actions_before_time_step2 ( const double &  time_step)
inlinevirtual

A virtual method that is run before the second time update.

For user-custom problems.

Parameters
time_stepThe time step size that is about to be taken

Definition at line 119 of file TwoD_TVDLF_Mesh.h.

119 {
120 // Empty by default. If you want to take actions before the time
121 // step then you need to inherit from this basic mesh and implement
122 // this method.
123 }

Referenced by update().

◆ boundary_diff()

void CppNoddy::TwoD_TVDLF_Mesh::boundary_diff ( elt_iter  e,
const int &  face_index,
DenseVector< double > &  diff 
) const
protected

Given an element iterator and the local coordinate this will return zero for any components specified as inflow boundary conditions in line with Levy & Tadmor (1997) and set the centre nodal value to the edge value.

Parameters
eElement iterator
face_indexThe index of the face that is on a boundary
diffA vector of derivatives

Definition at line 892 of file TwoD_TVDLF_Mesh.cpp.

892 {
893 DenseVector<double> se(2, 0.0);
894 // get the appropriate local coordinate for the mid point on the edge
895 switch(face_index) {
896 case 0:
897 // south
898 se[ 0 ] = 0.0;
899 se[ 1 ] = -1.0;
900 break;
901 case 1:
902 // east
903 se[ 0 ] = 1.0;
904 se[ 1 ] = 0.0;
905 break;
906 case 2:
907 // north
908 se[ 0 ] = 0.0;
909 se[ 1 ] = 1.0;
910 break;
911 case 3:
912 // west
913 se[ 0 ] = -1.0;
914 se[ 1 ] = 0.0;
915 break;
916 }
917 // for the edge value
918 DenseVector<double> Q(e -> get_Q(se));
919 // look at the user-defined BCs for any defined Dirichlet conditions
920 std::vector<bool> inflow = e -> p_system -> edge_values(face_index, e -> get_x(se), Q);
921 // the default is a zero slope in accordance with Levy & Tadmor (1997)
922 DenseVector<double> sigma_n(ORDER_OF_SYSTEM, 0.0);
923 // allow the user to override the zero slope
924 e -> p_system -> edge_slopes(face_index, e -> get_x(se), sigma_n);
925 // use edge values for inflow conditions
926 for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
927 if(inflow[ i ] == true) {
928 // set a zero slope as per Levy & Tadmor (1997)?
929 diff[ i ] = sigma_n[ i ];// ( Qe[ i ] - Qm[ i ] ) / delta;
930 //std::cout << face_index << " " << diff[ i ] << "\n";
931 }
932 }
933 }

References ORDER_OF_SYSTEM.

Referenced by east_diff(), EW_diff(), north_diff(), NS_diff(), south_diff(), and west_diff().

◆ calc_slopes()

void CppNoddy::TwoD_TVDLF_Mesh::calc_slopes ( vector_of_elts *  elt_vector)
protected

Use the appropriate limiter to approximate the slope in each element in the mesh.

The slopes will be sent down to the element objects and then onto the face objects.

Parameters
elt_vectorThe vector of elements to set the slope for

superbee

Definition at line 690 of file TwoD_TVDLF_Mesh.cpp.

690 {
691 set_boundary_Q(elt_vector);
692 elt_iter e(elt_vector -> begin());
693 DenseVector<double> slope_east(ORDER_OF_SYSTEM, 0.0);
694 DenseVector<double> slope_west(ORDER_OF_SYSTEM, 0.0);
695 DenseVector<double> slope_north(ORDER_OF_SYSTEM, 0.0);
696 DenseVector<double> slope_south(ORDER_OF_SYSTEM, 0.0);
697 switch(LIMITER) {
698 case 0:
699 // minmod
700 while(e != elt_vector -> end()) {
701 slope_east = east_diff(e);
702 slope_west = west_diff(e);
703 slope_south = south_diff(e);
704 slope_north = north_diff(e);
705 e -> set_slope_x(minmod(slope_east, slope_west));
706 e -> set_slope_y(minmod(slope_north, slope_south));
707 ++e;
708 }
709 break;
710 case 1:
711 // MC
712 while(e != elt_vector -> end()) {
713 slope_east = east_diff(e);
714 slope_west = west_diff(e);
715 DenseVector<double> slope_ew(EW_diff(e));
716 slope_south = south_diff(e);
717 slope_north = north_diff(e);
718 DenseVector<double> slope_ns(NS_diff(e));
719 const DenseVector<double> temp_x(minmod(slope_east * 2, slope_west * 2));
720 const DenseVector<double> slope_x(minmod(temp_x, slope_ew));
721 const DenseVector<double> temp_y(minmod(slope_north * 2, slope_south * 2));
722 const DenseVector<double> slope_y(minmod(temp_y, slope_ns));
723 e -> set_slope_x(slope_x);
724 e -> set_slope_y(slope_y);
725 ++e;
726 }
727 case 2:
728 /// superbee
729 while(e != elt_vector -> end()) {
730 slope_east = east_diff(e);
731 slope_west = west_diff(e);
732 slope_south = south_diff(e);
733 slope_north = north_diff(e);
734 const DenseVector<double> slope_x = maxmod(
735 minmod(slope_east, slope_west * 2.),
736 minmod(slope_east * 2., slope_west)
737 );
738 e -> set_slope_x(slope_x);
739 const DenseVector<double> slope_y = maxmod(
740 minmod(slope_north, slope_south * 2.),
741 minmod(slope_north * 2., slope_south)
742 );
743 e -> set_slope_y(slope_y);
744 ++e;
745 }
746 break;
747 default:
748 std::string problem;
749 problem = " The TwoD_TVDLF_Mesh object has an unrecognised 'LIMITER' identifier. \n";
750 throw ExceptionRuntime(problem);
751 }
752 }
CppNoddy::TwoD_TVDLF_Mesh::minmod
DenseVector< double > minmod(DenseVector< double > A, DenseVector< double > B) const
A vector version of the minmod operator.
Definition: TwoD_TVDLF_Mesh.cpp:874
CppNoddy::TwoD_TVDLF_Mesh::EW_diff
DenseVector< double > EW_diff(elt_iter e)
Compute a finite difference approximation of the derivative in the compass direction.
Definition: TwoD_TVDLF_Mesh.cpp:843
CppNoddy::TwoD_TVDLF_Mesh::west_diff
DenseVector< double > west_diff(elt_iter e) const
Compute a finite difference approximation of the derivative in the compass direction.
Definition: TwoD_TVDLF_Mesh.cpp:771
CppNoddy::TwoD_TVDLF_Mesh::east_diff
DenseVector< double > east_diff(elt_iter e) const
Compute a finite difference approximation of the derivative in the compass direction.
Definition: TwoD_TVDLF_Mesh.cpp:754
CppNoddy::TwoD_TVDLF_Mesh::south_diff
DenseVector< double > south_diff(elt_iter e) const
Compute a finite difference approximation of the derivative in the compass direction.
Definition: TwoD_TVDLF_Mesh.cpp:805
CppNoddy::TwoD_TVDLF_Mesh::set_boundary_Q
void set_boundary_Q(vector_of_elts *elts)
Loops over all boundary elements and sets the Q values in each one to be the value specified by the e...
Definition: TwoD_TVDLF_Mesh.h:319
CppNoddy::TwoD_TVDLF_Mesh::NS_diff
DenseVector< double > NS_diff(elt_iter e)
Compute a finite difference approximation of the derivative in the compass direction.
Definition: TwoD_TVDLF_Mesh.cpp:822
CppNoddy::TwoD_TVDLF_Mesh::maxmod
DenseVector< double > maxmod(DenseVector< double > A, DenseVector< double > B) const
A vector version of the maxmod operator.
Definition: TwoD_TVDLF_Mesh.cpp:883
CppNoddy::TwoD_TVDLF_Mesh::north_diff
DenseVector< double > north_diff(elt_iter e) const
Compute a finite difference approximation of the derivative in the compass direction.
Definition: TwoD_TVDLF_Mesh.cpp:788

References east_diff(), EW_diff(), LIMITER, maxmod(), minmod(), north_diff(), NS_diff(), ORDER_OF_SYSTEM, set_boundary_Q(), south_diff(), and west_diff().

Referenced by TwoD_TVDLF_Mesh(), and update().

◆ dump_data()

void CppNoddy::TwoD_TVDLF_Mesh::dump_data ( std::string  filename)

Dump the data over the nodal positions to a given filename.

Parameters
filenameThe filename to be generated

Definition at line 465 of file TwoD_TVDLF_Mesh.cpp.

465 {
466 std::ofstream dump;
467 dump.open(filename.c_str());
468 dump.precision(6);
469 elt_iter e;
470 DenseVector<double> s(2, 0.0);
471 {
472 e = BLACK_ELTS.begin();
473 while(e != BLACK_ELTS.end()) {
474 for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
475 dump << e -> get_Q(s)[i] << " ";
476 }
477 dump << e -> get_max_dt() << " ";
478 // for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
479 // dump << e -> get_slope_x()[i] << " ";
480 // }
481 // for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
482 // dump << e -> get_slope_y()[i] << " ";
483 // }
484 dump << "\n";
485 ++e;
486 }
487 }
488 dump.close();
489 }

References BLACK_ELTS, and ORDER_OF_SYSTEM.

◆ dump_gnu()

void CppNoddy::TwoD_TVDLF_Mesh::dump_gnu ( std::string  filename)

Dump the data to a given filename in a gnuplot format.

Parameters
filenameThe filename to be generated

Definition at line 405 of file TwoD_TVDLF_Mesh.cpp.

405 {
406 std::ofstream dump;
407 dump.open(filename.c_str());
408 dump.precision(6);
409 DenseVector<double> s(2, 0.0);
410 {
411 elt_iter e(BLACK_ELTS.begin());
412 while(e != BLACK_ELTS.end()) {
413 dump << e -> get_x(s)[0] << " " << e -> get_x(s)[1] << " ";
414 for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
415 dump << e -> get_Q(s)[i] << " ";
416 }
417 dump << e -> get_max_dt();
418 // for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
419 // dump << e -> get_slope_x()[i] << " ";
420 // }
421 // for(std::size_t i = 0; i < ORDER_OF_SYSTEM; ++i) {
422 // dump << e -> get_slope_y()[i] << " ";
423 // }
424 dump << "\n";
425 if(e -> face_is_external(1))
426 dump << "\n";
427 ++e;
428 }
429 }
430 dump.close();
431 }

References BLACK_ELTS, and ORDER_OF_SYSTEM.

Referenced by main().

◆ dump_nodes_x()

void CppNoddy::TwoD_TVDLF_Mesh::dump_nodes_x ( std::string  filename) const

Dump the x-nodal positions to a given filename.

This method only applies to regular structured meshes.

Parameters
filenameThe filename to be generated

Definition at line 433 of file TwoD_TVDLF_Mesh.cpp.

433 {
434 std::ofstream dump;
435 dump.open(filename.c_str());
436 dump.precision(6);
437 DenseVector<double> s(2, 0.0);
438 celt_iter e;
439 {
440 e = BLACK_ELTS.begin();
441 for(unsigned i = 0; i < NX - 1; ++i) {
442 dump << e -> get_x(s)[0] << "\n";
443 ++e;
444 }
445 }
446 dump.close();
447 }
CppNoddy::TwoD_TVDLF_Mesh::celt_iter
vector_of_elts::const_iterator celt_iter
Definition: TwoD_TVDLF_Mesh.h:23

References BLACK_ELTS, and NX.

◆ dump_nodes_y()

void CppNoddy::TwoD_TVDLF_Mesh::dump_nodes_y ( std::string  filename) const

Dump the y-nodal positions to a given filename.

This method only applies to regular structured meshes.

Parameters
filenameThe filename to be generated

Definition at line 449 of file TwoD_TVDLF_Mesh.cpp.

449 {
450 std::ofstream dump;
451 dump.open(filename.c_str());
452 dump.precision(6);
453 DenseVector<double> s(2, 0.0);
454 celt_iter e;
455 {
456 e = BLACK_ELTS.begin();
457 for(unsigned i = 0; i < NY - 1; ++i) {
458 dump << e -> get_x(s)[1] << "\n";
459 e += NX - 1;
460 }
461 }
462 dump.close();
463 }

References BLACK_ELTS, NX, and NY.

◆ east_diff()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::east_diff ( elt_iter  e) const
protected

Compute a finite difference approximation of the derivative in the compass direction.

Parameters
eThe iterator to the element that we want the derivative for
Returns
The spatial derivative of the vector system

Definition at line 754 of file TwoD_TVDLF_Mesh.cpp.

754 {
755 DenseVector<double> diff(ORDER_OF_SYSTEM, 0.0);
756 std::set< int > faces(e -> get_external_faces());
757 const int face_index(1);
758 if(e -> face_is_external(face_index)) {
759 // interior difference
760 diff = west_diff(e);
761 boundary_diff(e, face_index, diff);
762 } else {
763 // We're not on an edge, so we can difference
764 elt_iter j(e -> get_ptrs(face_index));
765 diff = (j -> get_Q_mid() - e -> get_Q_mid())
766 / (j -> get_x_mid()[0] - e -> get_x_mid()[0]);
767 }
768 return diff;
769 }
CppNoddy::TwoD_TVDLF_Mesh::boundary_diff
void boundary_diff(elt_iter e, const int &face_index, DenseVector< double > &diff) const
Given an element iterator and the local coordinate this will return zero for any components specified...
Definition: TwoD_TVDLF_Mesh.cpp:892

References boundary_diff(), ORDER_OF_SYSTEM, and west_diff().

Referenced by calc_slopes(), EW_diff(), and west_diff().

◆ EW_diff()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::EW_diff ( elt_iter  e)
protected

Compute a finite difference approximation of the derivative in the compass direction.

Parameters
eThe iterator to the element that we want the derivative for
Returns
The spatial derivative of the vector system

Definition at line 843 of file TwoD_TVDLF_Mesh.cpp.

843 {
844 DenseVector<double> diff(ORDER_OF_SYSTEM, 0.0);
845 std::set< int > faces(e -> get_external_faces());
846 if(e -> face_is_external(1)) {
847 // interior difference
848 diff = west_diff(e);
849 boundary_diff(e, 1, diff);
850 } else if(e -> face_is_external(3)) {
851 // interior difference
852 diff = east_diff(e);
853 boundary_diff(e, 3, diff);
854 } else {
855 // We're not on an edge, so we can difference
856 elt_iter jE(e -> get_ptrs(1));
857 elt_iter jW(e -> get_ptrs(3));
858 diff = (jE -> get_Q_mid() - jW -> get_Q_mid())
859 / (jE -> get_x_mid()[0] - jW -> get_x_mid()[0]);
860 }
861 return diff;
862 }

References boundary_diff(), east_diff(), ORDER_OF_SYSTEM, and west_diff().

Referenced by calc_slopes().

◆ get_elt_iter_from_elt_iter()

elt_iter CppNoddy::TwoD_TVDLF_Mesh::get_elt_iter_from_elt_iter ( elt_iter  e,
std::string  target_colour,
int  corner_index 
)
inlineprotected

Given an element in the INITIAL STRUCTURED MESH, this method will return an iterator to an element in the other mesh that overlaps the corner specified by corner_index.

Parameters
eAn element iterator to the source element
target_colourThe colour of the target mesh, ie. NOT the one that iterator 'e' is in
corner_indexThe index of the corner to be considered 0,1,2,3 for SW, SE, NE, NW.

Definition at line 139 of file TwoD_TVDLF_Mesh.h.

141 {
142 // for corners SW,SE,NE,NW:
143 // black (i,j) -> red (i,j),(i+1,j),(i+1,j+1),(i,j+1)
144 // => j*Bx+i -> j*Rx+i etc
145 // red (i,j) -> black (i-1,j-1),(i,j-1),(i,j),(i-1,j)
146 // => j*Rx+i -> (j-1)*Bx+i-1 etc
147 //dimensions of the Black (Bx times By) & Red (Rx times Ry) meshes
148 //
149 vector_of_elts* target_elts(NULL);
150 int k(0);
151 int Bx(NX - 1);
152 int By(NY - 1);
153 int Rx(NX);
154 int Ry(NY);
155 //
156 if(target_colour == "red") {
157 int offset(e - BLACK_ELTS.begin());
158 target_elts = &RED_ELTS;
159 int i, j, l, m;
160 //std::cout << " offset = " << offset << "\n";
161 i = offset % Bx;
162 j = (offset - i) / Bx;
163 switch(corner_index) {
164 case 0:
165 l = i;
166 m = j;
167 break;
168 case 1:
169 l = i + 1;
170 m = j;
171 break;
172 case 2:
173 l = i + 1;
174 m = j + 1;
175 break;
176 case 3:
177 l = i;
178 m = j + 1;
179 break;
180 default:
181 l = 0;
182 m = 0;
183 assert(false);
184 }
185 // dont return outside the mesh
186 //std::cout << source_colour << " " << i << "," << j << " -> " << l << "," << m <<"\n";
187 l = std::min(l, Rx);
188 l = std::max(l, 0);
189 m = std::min(m, Ry);
190 m = std::max(m, 0);
191 k = m * Rx + l;
192 //std::cout << source_colour << " " << i << "," << j << " -> " << l << "," << m <<"\n";
193 //std::cout << "--\n";
194 } else if(target_colour == "black") {
195 int offset(e - RED_ELTS.begin());
196 target_elts = &BLACK_ELTS;
197 int i, j, l, m;
198 i = offset % Rx;
199 j = (offset - i) / Rx;
200 switch(corner_index) {
201 case 0:
202 l = i - 1;
203 m = j - 1;
204 break;
205 case 1:
206 l = i;
207 m = j - 1;
208 break;
209 case 2:
210 l = i;
211 m = j;
212 break;
213 case 3:
214 l = i - 1;
215 m = j;
216 break;
217 default:
218 l = 0;
219 m = 0;
220 assert(false);
221 }
222 // dont return outside the mesh
223 l = std::min(l, Bx);
224 l = std::max(l, 0);
225 m = std::min(m, By);
226 m = std::max(m, 0);
227 k = m * Bx + l;
228 } else {
229 // throw
230 }
231 return target_elts -> begin() + k;
232 }
CppNoddy::TwoD_TVDLF_Mesh::vector_of_elts
std::vector< TwoD_TVDLF_Elt > vector_of_elts
iterators for the vector of elements
Definition: TwoD_TVDLF_Mesh.h:22

References BLACK_ELTS, m, NX, NY, and RED_ELTS.

Referenced by TwoD_TVDLF_Mesh().

◆ get_elt_iter_from_x()

std::vector< TwoD_TVDLF_Elt >::iterator CppNoddy::TwoD_TVDLF_Mesh::get_elt_iter_from_x ( const DenseVector< double > &  x,
std::string  mesh_colour = "black" 
)

Given a global coordinate, return a pointer to the elt that contains that point.

Because the elts are stored in a vector, and are logically in a rectangular grid, we first hunt for the column, then the appropriate row.

Parameters
xA vector global coordinate
mesh_colourA string identifier of the mesh to be searched
Returns
A pointer to the element containing the global coordinate

Definition at line 618 of file TwoD_TVDLF_Mesh.cpp.

618 {
619 elt_iter found_elt;
620 vector_of_elts* elts(get_elts_from_colour(mesh_colour));
621 std::size_t Mx(get_number_elts_in_x(mesh_colour));
622
623 elt_iter e = elts -> begin();
624 while(e != elts -> end()) {
625 // local coordinates of the NE and SW corners
626 DenseVector<double> ne(2, 1.0);
627 DenseVector<double> sw(-ne);
628 double west = e -> get_x(sw)[ 0 ];
629 double east = e -> get_x(ne)[ 0 ];
630 if((x[ 0 ] >= west) && (x[ 0 ] <= east)) {
631 break;
632 }
633 ++e;
634 }
635
636 while(true) {
637 // local coordinates of the NE and SW corners
638 DenseVector<double> ne(2, 1.0);
639 DenseVector<double> sw(-ne);
640 double south = e -> get_x(sw)[ 1 ];
641 double north = e -> get_x(ne)[ 1 ];
642 if((x[ 1 ] >= south) && (x[ 1 ] <= north)) {
643 break;
644 }
645 if(e + Mx < elts -> end()) {
646 e += Mx;
647 } else {
648 // if we are here, then we're looking for a point outside the mesh
649 std::string problem;
650 problem = "The TwoD_TVDLF_Mesh::get_elt_iter_from_x method has been called\n";
651 problem += "for a position that is outside the boundaries of the mesh.\n";
652 throw ExceptionRuntime(problem);
653 }
654 }
655 // if we get to here, then we've found the elt.
656 return e;
657 }
CppNoddy::TwoD_TVDLF_Mesh::get_elts_from_colour
vector_of_elts * get_elts_from_colour(std::string mesh_colour)
Given a choice of black or red mesh, return a pointer to the appropriate vector of elements.
Definition: TwoD_TVDLF_Mesh.cpp:660
CppNoddy::TwoD_TVDLF_Mesh::get_number_elts_in_x
std::size_t get_number_elts_in_x(std::string mesh_colour)
Given a choice of black or red mesh, return the number of elements in the x-direction.
Definition: TwoD_TVDLF_Mesh.cpp:675

References get_elts_from_colour(), and get_number_elts_in_x().

Referenced by get_point_values().

◆ get_elts_from_colour()

std::vector< TwoD_TVDLF_Elt > * CppNoddy::TwoD_TVDLF_Mesh::get_elts_from_colour ( std::string  mesh_colour)
protected

Given a choice of black or red mesh, return a pointer to the appropriate vector of elements.

Parameters
mesh_colourA mesh identifier (black or red)
Returns
A pointer to a std::vector of elements

Definition at line 660 of file TwoD_TVDLF_Mesh.cpp.

660 {
661 vector_of_elts* elts;
662 if(mesh_colour == "black") {
663 elts = &BLACK_ELTS;
664 } else if(mesh_colour == "red") {
665 elts = &RED_ELTS;
666 } else {
667 std::string problem;
668 problem = " The TwoD_TVDLF_Mesh object has been passed an unrecognised mesh identifier. \n";
669 problem += " valid options are 'black' or 'red'.\n";
670 throw ExceptionRuntime(problem);
671 }
672 return elts;
673 }

References BLACK_ELTS, and RED_ELTS.

Referenced by get_elt_iter_from_x(), and integrate().

◆ get_number_elts_in_x()

std::size_t CppNoddy::TwoD_TVDLF_Mesh::get_number_elts_in_x ( std::string  mesh_colour)
protected

Given a choice of black or red mesh, return the number of elements in the x-direction.

Parameters
mesh_colourA mesh identifier (black or red)
Returns
The number of elts in the x-direction

Definition at line 675 of file TwoD_TVDLF_Mesh.cpp.

675 {
676 std::size_t N;
677 if(mesh_colour == "black") {
678 N = NX - 1;
679 } else if(mesh_colour == "red") {
680 N = NX;
681 } else {
682 std::string problem;
683 problem = " The TwoD_TVDLF_Mesh object has been passed an unrecognised mesh identifier. \n";
684 problem += " valid options are 'black' or 'red'.\n";
685 throw ExceptionRuntime(problem);
686 }
687 return N;
688 }

References NX.

Referenced by get_elt_iter_from_x().

◆ get_point_values()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::get_point_values ( const DenseVector< double > &  x)

Get the vector of unknowns at a point in the 2D mesh.

Parameters
xThe global coordinate at which to return the values
Returns
The vector of unknowns

Definition at line 491 of file TwoD_TVDLF_Mesh.cpp.

491 {
492 elt_iter e(get_elt_iter_from_x(x));
493 DenseVector<double> s(e-> get_s(x));
494 return e -> get_Q(s);
495 }
CppNoddy::TwoD_TVDLF_Mesh::get_elt_iter_from_x
elt_iter get_elt_iter_from_x(const DenseVector< double > &x, std::string mesh_colour="black")
Given a global coordinate, return a pointer to the elt that contains that point.
Definition: TwoD_TVDLF_Mesh.cpp:618

References get_elt_iter_from_x().

Referenced by main().

◆ get_time()

const double & CppNoddy::TwoD_TVDLF_Mesh::get_time ( ) const

Get a const reference to the time value for the current mesh.

Returns
time The time level at which the data in the mesh applies.

Definition at line 601 of file TwoD_TVDLF_Mesh.cpp.

601 {
602 return MESH_TIME;
603 }

References MESH_TIME.

Referenced by main().

◆ integrate()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::integrate ( std::string  mesh_colour = "black")

Integrate the concentration values across the entire mesh.

Parameters
mesh_colourThe identifier of which mesh to integrate over
Returns
The vector value of the integral

Definition at line 605 of file TwoD_TVDLF_Mesh.cpp.

605 {
606 vector_of_elts* elts(get_elts_from_colour(mesh_colour));
607 elt_iter e = elts -> begin();
608 DenseVector<double> I(ORDER_OF_SYSTEM, 0.0);
609 while(e != elts -> end()) {
610 const DenseVector<double> sw(2, -1.0);
611 const DenseVector<double> ne(2, 1.0);
612 I += e -> get_int_Q(sw, ne);
613 ++e;
614 }
615 return I;
616 }

References get_elts_from_colour(), and ORDER_OF_SYSTEM.

◆ maxmod()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::maxmod ( DenseVector< double >  A,
DenseVector< double >  B 
) const
protected

A vector version of the maxmod operator.

Parameters
AA vector to compare
BA vector to compare
Returns
A component-wise maxmod vector, i.e., each component of the returned vector is the maxmod of the components of A & B.

Definition at line 883 of file TwoD_TVDLF_Mesh.cpp.

883 {
884 DenseVector<double> MM;
885 for(std::size_t i = 0; i < A.size(); ++i) {
886 MM.push_back(0.5 * (sgn(A[i]) + sgn(B[i]))
887 * std::max(std::abs(A[i]), std::abs(B[i])));
888 }
889 return MM;
890 }
CppNoddy::TwoD_TVDLF_Mesh::sgn
int sgn(double a) const
Sign of a double.
Definition: TwoD_TVDLF_Mesh.cpp:864
CppNoddy::Example::A
double A(1.0)
initial hump amplitude

References CppNoddy::DenseVector< _Type >::push_back(), and sgn().

Referenced by calc_slopes().

◆ minmod()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::minmod ( DenseVector< double >  A,
DenseVector< double >  B 
) const
protected

A vector version of the minmod operator.

Parameters
AA vector to compare
BA vector to compare
Returns
A component-wise minmod vector, i.e., each component of the returned vector is the minmod of the components of A & B.

Definition at line 874 of file TwoD_TVDLF_Mesh.cpp.

874 {
875 DenseVector<double> MM;
876 for(std::size_t i = 0; i < A.size(); ++i) {
877 MM.push_back(0.5 * (sgn(A[i]) + sgn(B[i]))
878 * std::min(std::abs(A[i]), std::abs(B[i])));
879 }
880 return MM;
881 }

References CppNoddy::DenseVector< _Type >::push_back(), and sgn().

Referenced by calc_slopes().

◆ north_diff()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::north_diff ( elt_iter  e) const
protected

Compute a finite difference approximation of the derivative in the compass direction.

Parameters
eThe iterator to the element that we want the derivative for
Returns
The spatial derivative of the vector system

Definition at line 788 of file TwoD_TVDLF_Mesh.cpp.

788 {
789 DenseVector<double> diff(ORDER_OF_SYSTEM, 0.0);
790 std::set< int > faces(e -> get_external_faces());
791 const int face_index(2);
792 if(e -> face_is_external(face_index)) {
793 // interior difference
794 diff = south_diff(e);
795 boundary_diff(e, face_index, diff);
796 } else {
797 // We're not on an edge, so we can difference
798 elt_iter j(e -> get_ptrs(face_index));
799 diff = (j -> get_Q_mid() - e -> get_Q_mid())
800 / (j -> get_x_mid()[1] - e -> get_x_mid()[1]);
801 }
802 return diff;
803 }

References boundary_diff(), ORDER_OF_SYSTEM, and south_diff().

Referenced by calc_slopes(), NS_diff(), and south_diff().

◆ NS_diff()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::NS_diff ( elt_iter  e)
protected

Compute a finite difference approximation of the derivative in the compass direction.

Parameters
eThe iterator to the element that we want the derivative for
Returns
The spatial derivative of the vector system

Definition at line 822 of file TwoD_TVDLF_Mesh.cpp.

822 {
823 DenseVector<double> diff(ORDER_OF_SYSTEM, 0.0);
824 std::set< int > faces(e -> get_external_faces());
825 if(e -> face_is_external(0)) {
826 // interior difference
827 diff = north_diff(e);
828 boundary_diff(e, 0, diff);
829 } else if(e -> face_is_external(2)) {
830 // interior difference
831 diff = south_diff(e);
832 boundary_diff(e, 2, diff);
833 } else {
834 // We're not on an edge, so we can difference
835 elt_iter jS(e -> get_ptrs(0));
836 elt_iter jN(e -> get_ptrs(2));
837 diff = (jN -> get_Q_mid() - jS -> get_Q_mid())
838 / (jN -> get_x_mid()[1] - jS -> get_x_mid()[1]);
839 }
840 return diff;
841 }

References boundary_diff(), north_diff(), ORDER_OF_SYSTEM, and south_diff().

Referenced by calc_slopes().

◆ set_boundary_Q()

void CppNoddy::TwoD_TVDLF_Mesh::set_boundary_Q ( vector_of_elts *  elts)
inlineprotected

Loops over all boundary elements and sets the Q values in each one to be the value specified by the edge_values method.

Parameters
eltsA vector of elements in the mesh

Definition at line 319 of file TwoD_TVDLF_Mesh.h.

319 {
320 // loop over all elts
321 elt_iter e(elts -> begin());
322 while(e != elts -> end()) {
323 // find the external elts that are on the boundary
324 if(e -> get_external_flag()) {
325 std::set<int> faces(e -> get_external_faces());
326 std::set<int>::iterator i(faces.begin());
327 // loop through all external faces (there are 2 on corner elts)
328 while(i != faces.end()) {
329 int face(*i);
330 // local coordinate at which to impose the BC
331 DenseVector<double> s(2, 0.0);
332 switch(face) {
333 case 0:
334 s[ 0 ] = 0.0;
335 s[ 1 ] = -1.0;
336 break;
337 case 1:
338 s[ 0 ] = 1.0;
339 s[ 1 ] = 0.0;
340 break;
341 case 2:
342 s[ 0 ] = 0.0;
343 s[ 1 ] = 1.0;
344 break;
345 case 3:
346 s[ 0 ] = -1.0;
347 s[ 1 ] = 0.0;
348 break;
349 }
350 // get the current edge data
351 DenseVector<double> Qe(e -> get_Q(s));
352 //// initialise the normal slope vector
353 // DenseVector<double> sigma_n( ORDER_OF_SYSTEM, 0.0 );
354 // get the user specified edge values
355 std::vector<bool> inflow = e -> p_system -> edge_values(face, e -> get_x(s), Qe); //, sigma_n );
356 // get the mid-element nodal data
357 DenseVector<double> Qm(e -> get_Q_mid());
358 for(std::size_t j = 0; j < ORDER_OF_SYSTEM; ++j) {
359 if(inflow[ j ] == true) {
360 // only change components that are inflow
361 Qm[ j ] = Qe[ j ];
362 }
363 }
364 e -> set_Q_mid(Qm);
365 ++i;
366 }
367 }
368 ++e;
369 }
370 }

References ORDER_OF_SYSTEM.

Referenced by calc_slopes().

◆ set_limiter()

void CppNoddy::TwoD_TVDLF_Mesh::set_limiter ( const unsigned &  id)

Set the limiter type to be applied in the slope values.

Parameters
idThe identifier of the limiter.

Definition at line 497 of file TwoD_TVDLF_Mesh.cpp.

497 {
498 LIMITER = id;
499 }

References LIMITER.

Referenced by main().

◆ sgn()

int CppNoddy::TwoD_TVDLF_Mesh::sgn ( double  a) const
protected

Sign of a double.

Parameters
aThe value to return the sign of
Returns
The sign of the value

Definition at line 864 of file TwoD_TVDLF_Mesh.cpp.

864 {
865 if(a > 0.0) {
866 return 1;
867 } else if(a < 0.0) {
868 return -1;
869 } else {
870 return 0;
871 }
872 }

Referenced by maxmod(), and minmod().

◆ south_diff()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::south_diff ( elt_iter  e) const
protected

Compute a finite difference approximation of the derivative in the compass direction.

Parameters
eThe iterator to the element that we want the derivative for
Returns
The spatial derivative of the vector system

Definition at line 805 of file TwoD_TVDLF_Mesh.cpp.

805 {
806 DenseVector<double> diff(ORDER_OF_SYSTEM, 0.0);
807 std::set< int > faces(e -> get_external_faces());
808 const int face_index(0);
809 if(e -> face_is_external(face_index)) {
810 // interior difference
811 diff = north_diff(e);
812 boundary_diff(e, face_index, diff);
813 } else {
814 // We're not on an edge, so we can difference
815 elt_iter j(e -> get_ptrs(face_index));
816 diff = (e -> get_Q_mid() - j -> get_Q_mid())
817 / (e -> get_x_mid()[1] - j -> get_x_mid()[1]);
818 }
819 return diff;
820 }

References boundary_diff(), north_diff(), and ORDER_OF_SYSTEM.

Referenced by calc_slopes(), north_diff(), and NS_diff().

◆ update()

double CppNoddy::TwoD_TVDLF_Mesh::update ( const double &  CFL,
const double &  max_dt = std::numeric_limits<long double>::max() 
)

Update the mesh object.

A time step is chosen based upon the CFL constraint.

Parameters
CFLThe CFL constraint for the timestep (e.g. 0.49)
max_dtA maximum timestep to be taken irrespective of the CFL value
Returns
The total timestep that was taken

Definition at line 501 of file TwoD_TVDLF_Mesh.cpp.

501 {
502 // integrate the black mesh data onto the red mesh
503 std::cout << "[DEBUG] in update with MESH_TIME = " << MESH_TIME << "\n";
504 std::cout << "[DEBUG] in update with max_dt = " << max_dt << "\n";
505
506 {
507 elt_iter e = RED_ELTS.begin();
508 while(e != RED_ELTS.end()) {
509 e -> set_Q_mid(e -> contributed_Q() / e -> get_dA());
510 ++e;
511 }
512 }
513 calc_slopes(&RED_ELTS);
514
515 // determine the first time step from the CFL constraint
516 double first_dt;
517 {
518 elt_iter e = BLACK_ELTS.begin();
519 first_dt = e -> get_max_dt();
520 ++e;
521 while(e != BLACK_ELTS.end()) {
522 first_dt = std::min(first_dt, e -> get_max_dt());
523 ++e;
524 }
525 std::cout << "[DEBUG] in update with first_dt = " << first_dt << "\n";
526 first_dt *= CFL;
527 std::cout << "[DEBUG] in update with first_dt*CFL = " << first_dt << "\n";
528 }
529 if(first_dt > max_dt / 2) {
530 first_dt = max_dt / 2;
531 }
532 std::cout << "[DEBUG] in update with final first_dt = " << first_dt << "\n";
533
534 actions_before_time_step1(first_dt);
535 // do the time step
536 {
537 elt_iter e = RED_ELTS.begin();
538 while(e != RED_ELTS.end()) {
539 // add to the corrections
540 e -> add_flux_contributions(first_dt);
541 ++e;
542 }
543 }
544 calc_slopes(&RED_ELTS);
545 MESH_TIME += first_dt;
546 std::cout << "[DEBUG] in update with MESH_TIME = " << MESH_TIME << "\n";
547
548 // integrate the red mesh data onto the black mesh
549 {
550 elt_iter e = BLACK_ELTS.begin();
551 while(e != BLACK_ELTS.end()) {
552 e -> set_Q_mid(e -> contributed_Q() / e -> get_dA());
553 ++e;
554 }
555 }
556 calc_slopes(&BLACK_ELTS);
557
558 // determine the first time step from the CFL constraint
559 double second_dt;
560 {
561 elt_iter e = RED_ELTS.begin();
562 second_dt = e -> get_max_dt();
563 ++e;
564 while(e != RED_ELTS.end()) {
565 second_dt = std::min(second_dt, e -> get_max_dt());
566 ++e;
567 }
568 std::cout << "[DEBUG] in update with second_dt = " << second_dt << "\n";
569 second_dt *= CFL;
570 std::cout << "[DEBUG] in update with second_dt*CFL = " << second_dt << "\n";
571 }
572 if(first_dt + second_dt > max_dt) {
573 second_dt = max_dt - first_dt;
574 }
575 std::cout << "[DEBUG] in update with final second_dt = " << second_dt << "\n";
576
577 actions_before_time_step2(second_dt);
578 // do the time step
579 {
580 elt_iter e = BLACK_ELTS.begin();
581 while(e != BLACK_ELTS.end()) {
582 // add to the corrections
583 e -> add_flux_contributions(second_dt);
584 ++e;
585 }
586 }
587 calc_slopes(&BLACK_ELTS);
588 MESH_TIME += second_dt;
589 std::cout << "[DEBUG] in update with MESH_TIME = " << MESH_TIME << "\n";
590
591 return first_dt + second_dt;
592 }
CppNoddy::TwoD_TVDLF_Mesh::actions_before_time_step1
virtual void actions_before_time_step1(const double &time_step)
A virtual method that is run before the first time update.
Definition: TwoD_TVDLF_Mesh.h:110
CppNoddy::TwoD_TVDLF_Mesh::actions_before_time_step2
virtual void actions_before_time_step2(const double &time_step)
A virtual method that is run before the second time update.
Definition: TwoD_TVDLF_Mesh.h:119

References actions_before_time_step1(), actions_before_time_step2(), BLACK_ELTS, calc_slopes(), MESH_TIME, and RED_ELTS.

Referenced by main(), and update_to().

◆ update_to()

void CppNoddy::TwoD_TVDLF_Mesh::update_to ( const double &  CFL,
const double &  t_end 
)

Update the mesh object to a set time level.

A time step is chosen based upon the CFL constraint.

Parameters
CFLThe CFL constraint for the timestep (e.g. 0.49)
t_endThe target end time.

Definition at line 594 of file TwoD_TVDLF_Mesh.cpp.

594 {
595 do {
596 std::cout << " computing to " << t_end << "\n";
597 update(CFL, std::abs(t_end - MESH_TIME));
598 } while(MESH_TIME < t_end);
599 }
CppNoddy::TwoD_TVDLF_Mesh::update
double update(const double &CFL, const double &max_dt=std::numeric_limits< long double >::max())
Update the mesh object.
Definition: TwoD_TVDLF_Mesh.cpp:501

References MESH_TIME, and update().

◆ west_diff()

DenseVector< double > CppNoddy::TwoD_TVDLF_Mesh::west_diff ( elt_iter  e) const
protected

Compute a finite difference approximation of the derivative in the compass direction.

Parameters
eThe iterator to the element that we want the derivative for
Returns
The spatial derivative of the vector system

Definition at line 771 of file TwoD_TVDLF_Mesh.cpp.

771 {
772 DenseVector<double> diff(ORDER_OF_SYSTEM, 0.0);
773 std::set< int > faces(e -> get_external_faces());
774 const int face_index(3);
775 if(e -> face_is_external(face_index)) {
776 // interior difference
777 diff = east_diff(e);
778 boundary_diff(e, face_index, diff);
779 } else {
780 // We're not on an edge, so we can difference
781 elt_iter j(e -> get_ptrs(face_index));
782 diff = (e -> get_Q_mid() - j -> get_Q_mid())
783 / (e -> get_x_mid()[0] - j -> get_x_mid()[0]);
784 }
785 return diff;
786 }

References boundary_diff(), east_diff(), and ORDER_OF_SYSTEM.

Referenced by calc_slopes(), east_diff(), and EW_diff().

Member Data Documentation

◆ BLACK_ELTS

vector_of_elts CppNoddy::TwoD_TVDLF_Mesh::BLACK_ELTS

An STL vector of linear elements – the black mesh.

Definition at line 126 of file TwoD_TVDLF_Mesh.h.

Referenced by dump_data(), dump_gnu(), dump_nodes_x(), dump_nodes_y(), get_elt_iter_from_elt_iter(), get_elts_from_colour(), TwoD_TVDLF_Mesh(), and update().

◆ LIMITER

unsigned CppNoddy::TwoD_TVDLF_Mesh::LIMITER
protected

Slope limiter method.

Definition at line 373 of file TwoD_TVDLF_Mesh.h.

Referenced by calc_slopes(), set_limiter(), and TwoD_TVDLF_Mesh().

◆ MESH_TIME

double CppNoddy::TwoD_TVDLF_Mesh::MESH_TIME
protected

the time level of the mesh

Definition at line 379 of file TwoD_TVDLF_Mesh.h.

Referenced by get_time(), TwoD_TVDLF_Mesh(), update(), and update_to().

◆ NX

std::size_t CppNoddy::TwoD_TVDLF_Mesh::NX
protected

number of faces in the x & y directions

Definition at line 377 of file TwoD_TVDLF_Mesh.h.

Referenced by dump_nodes_x(), dump_nodes_y(), get_elt_iter_from_elt_iter(), get_number_elts_in_x(), and TwoD_TVDLF_Mesh().

◆ NY

std::size_t CppNoddy::TwoD_TVDLF_Mesh::NY
protected

Definition at line 377 of file TwoD_TVDLF_Mesh.h.

Referenced by dump_nodes_y(), get_elt_iter_from_elt_iter(), and TwoD_TVDLF_Mesh().

◆ ORDER_OF_SYSTEM

std::size_t CppNoddy::TwoD_TVDLF_Mesh::ORDER_OF_SYSTEM
protected

order of the conservative system

Definition at line 375 of file TwoD_TVDLF_Mesh.h.

Referenced by boundary_diff(), calc_slopes(), dump_data(), dump_gnu(), east_diff(), EW_diff(), integrate(), north_diff(), NS_diff(), set_boundary_Q(), south_diff(), TwoD_TVDLF_Mesh(), and west_diff().

◆ p_Q_INIT

fn_ptr CppNoddy::TwoD_TVDLF_Mesh::p_Q_INIT
protected

function pointer to a funnction that defines the initial distribution

Definition at line 381 of file TwoD_TVDLF_Mesh.h.

Referenced by TwoD_TVDLF_Mesh().

◆ RED_ELTS

vector_of_elts CppNoddy::TwoD_TVDLF_Mesh::RED_ELTS

An STL vector of linear elements – the red mesh.

Definition at line 128 of file TwoD_TVDLF_Mesh.h.

Referenced by get_elt_iter_from_elt_iter(), get_elts_from_colour(), TwoD_TVDLF_Mesh(), and update().

The documentation for this class was generated from the following files:
Generated on Tue Jan 1 1980 00:00:00 for CppNoddy by   doxygen 1.9.5

© 2012

R.E. Hewitt