VERB_code_2.2  2
 All Classes Namespaces Files Functions Variables Typedefs Macros Pages
BoundaryConditions.cpp
Go to the documentation of this file.
1 /**
2  * \file BoundaryConditions.cpp
3  * Boundary condition class.
4  *
5  * \author Developed under supervision of the PI Yuri Shprits
6  */
7 
8 #include "BoundaryConditions.h"
9 #include <math.h>
10 
11 // Function 'strncasecmp' called '_strnicmp' in windows
12 #if defined(_WIN32) || defined(_WIN64)
13 #define strncasecmp _strnicmp
14 #define strcasecmp _stricmp
15 #endif
16 
17 
18 
19 /**
20  * Making boundary conditions.
21  *
22  * Fill-out 2D array of boundary conditions. (in general, 3D grid case) And save boundary condition type for future use.
23  *
24  * \param psd2DSlice - 2D array, slice of PSD on a boundary
25  * \param gridElement1 - one side of the grid
26  * \param gridElement2 - second side of the grid
27  */
29  int ix, iy;
30 
31  if (!arr.initialized) {
32  // Allocate memory
33  arr.AllocateMemory(psd2DSlice.size_x, psd2DSlice.size_y);
34  }
35 
36  // store type in a class variable, for the ease of future use
37  this->type = parameters.type;
38  if (parameters.type == "BCT_CONSTANT_VALUE") {
39  // Saving boundary value inside the class.
40  // If boundary is on constant value of the function - saving this type and the value.
41  this->calculationType = "BCT_CONSTANT_VALUE";
42 
43  // If value is zero, we don't put zero to PSD, but put some small value
44  //(*this) = VF::max(value, VC::zero_f);
45  arr = VF::max(parameters.value, VC::zero_f);
46 
47  } else if (parameters.type == "BCT_CONSTANT_PSD") {
48  // If boundary is on constant PSD value - saving this type and we need to store the initial values of PSD later, after initialization of PSD.
49  this->calculationType = "BCT_CONSTANT_VALUE";
50 
51  // So, this line store psd slice as a boundry, and change all zeros to some small value
52  for (ix = 0; ix < arr.size_x; ix++) {
53  for (iy = 0; iy < arr.size_y; iy++) {
54  arr[ix][iy] = VF::max(psd2DSlice[ix][iy], VC::zero_f);
55  }
56  }
57  } else if (parameters.type == "BCT_CONSTANT_DERIVATIVE") {
58  // If boundary is on constant derivative of the function - saving this type and the value of the derivative (usually - zero derivative).
59  this->calculationType = "BCT_CONSTANT_DERIVATIVE";
60 
61  // boundary condition is a constant derivative - so, storing derivative value for boundary conditions for each boundary point
62  // storing it inside the parent class Matrix.
63  // this - is a pointer to current class. '=' - operator.
64  arr = parameters.value;
65  } else if (parameters.type == "BCT_CONSTANT_SECOND_DERIVATIVE") {
66  // If boundary is on constant second derivative of the function
67  this->calculationType = "BCT_CONSTANT_SECOND_DERIVATIVE";
68 
69  // boundary condition is a constant derivative - so, storing derivative value for boundary conditions for each boundary point
70  // storing it inside the parent class Matrix.
71  // this - is a pointer to current class. '=' - operator.
72  arr = parameters.value;
73  } else if (parameters.type == "BCT_FILE" || parameters.type == "BCT_FILE_GRID" || parameters.type == "BCT_MULTIPLE_FILES") {
74  // If boundary are loaded from file... see BoundaryCondition::LoadBoundaryCondition
75  this->calculationType = "BCT_CONSTANT_VALUE";
76  // this->filename = parameters.filename; // don't need this i guess... but it does not do any bad.
77  } else {
78  throw error_msg("BCT_ERROR", "Unknown boundary condition type: %d", parameters.type.c_str());
79  }
80 }
81 
82 
83 /**
84  * Load boundary conditions from file.
85  *
86  * Here we actually fill-out 2D array of boundary conditions. (in general, 3D grid, case)
87  *
88  * \param Matrix2D<double> psd2DSlice - 2D array, slice of PSD on a boundary
89  * \param gridElement1 - one side of the grid
90  * \param gridElement2 - second side of the grid
91  */
93  int ix, iy;
94  double err = 1.e-3;
95  string inBuf;
96 
97  if (!arr.initialized) {
98  // Allocating memory
99  arr.AllocateMemory(gridElement1.size_x, gridElement1.size_y);
100  }
101 
102  Output::echo(0, "Loading boundary conditions from files... ");
103  // store type in a class variable, for the ease of future use
104  this->type = parameters.type;
105  //if (initialType == "BCT_FILE" || initialType == "BCT_FILE_GRID") {
106  if (parameters.type == "BCT_FILE" || parameters.type == "BCT_FILE_GRID") {
107  this->calculationType = "BCT_CONSTANT_VALUE";
108  // this->readFromFile(filename);
109  // making stream
110  {
111  ifstream input(parameters.filename.c_str());
112  // checking if file was opened
113  if (input == NULL) {
114  throw error_msg("BC_FILE_OPEN_ERR", "Boundary conditions file %s not found.", parameters.filename.c_str());
115  }
116  // Skip header!
117  input >> inBuf;
118  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) input >> inBuf;
119  input.ignore(9999, '\n'); // read to the end of the line with 'zone'
120  double loadedTime, loadedAx1, loadedAx2, loadedValue;
121  // size_x - time steps, size_y and size_z - two dimensions of PSD slice
122  for (ix = 0; ix < arr.size_x; ix++) {
123  for (iy = 0; iy < arr.size_y; iy++) {
124  if (parameters.type == "BCT_FILE_GRID") input >> loadedAx1 >> loadedAx2;
125  input >> loadedValue;
126  // check if grid is the same
127  if (parameters.type == "BCT_FILE_GRID" && (fabs(loadedAx1 - gridElement1[ix][iy]) > err || fabs(loadedAx2 - gridElement2[ix][iy]) > err)) {
128  if (input.eof()) {
129  throw error_msg("BC_LOAD_GRID_ERR", "Loading %s: not enough data points in the boundary condition file.\n", parameters.filename.c_str());
130  } else {
131  throw error_msg("BC_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: ax1=%e, ax2=%e\nGrid: ax1=%e, ax2=%e\n(i1=%d, i2=%d)\n", parameters.filename.c_str(), loadedAx1, loadedAx2, gridElement1[ix][iy], gridElement2[ix][iy], ix, iy);
132  }
133  } else {
134  arr[ix][iy] = loadedValue;
135  }
136  // skip till the end of the line
137  input.ignore(9999, '\n');
138  }
139  }
140  input.close();
141  }
142  arr.writeToFile("./Debug/" + arr.name + "_BC.dat");
143  Output::echo(0, "done.\n");
144  } else if (parameters.type == "BCT_MULTIPLE_FILES") {
145  // In that case we need to change boundary conditions (BC) with time.
146  // Now we load the list of boundary conditions (files with boundary values) and corresponding times, when we need to apply the BC
147  this->calculationType = "BCT_CONSTANT_VALUE";
148 
149  Output::echo(0, "Loading a list of boundary conditions... ");
150  ifstream input(parameters.filename.c_str());
151  // checking if file was opened
152  if (input == NULL) {
153  throw error_msg("BC_FILE_OPEN_ERR", "Boundary conditions file %s not found.", parameters.filename.c_str());
154  }
155 
156  // Skip header of the file
157  input >> inBuf;
158  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) input >> inBuf;
159  input.ignore(9999, '\n'); // read to the end of the line with 'zone'
160 
161  // reading the file
162  double time_tmp;
163  string filename_tmp;
164  while (!input.eof()) {
165  // reading time values and function values from the file into vectors
166  input >> time_tmp >> filename_tmp;
167  BC_change_time.push_back(time_tmp);
168  BC_filename.push_back(filename_tmp);
169  }
170 
171  // Load the first on
172  ifstream input2(BC_filename[0].c_str());
173  // checking if file was opened
174  if (input2 == NULL) {
175  throw error_msg("BC_FILE_OPEN_ERR", "Boundary conditions file %s not found.", parameters.filename.c_str());
176  }
177  Output::echo(3, "Updating a boundary condition: %s", BC_filename[0].c_str());
178  // Skip header of the file
179  string inBuf;
180  input2 >> inBuf;
181  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input2.eof() ) input2 >> inBuf;
182  input2.ignore(9999, '\n'); // read to the end of the line with 'zone'
183  double tmp;
184  for (ix = 0; ix < arr.size_x; ix++) {
185  for (iy = 0; iy < arr.size_y; iy++) {
186  input2 >> arr[ix][iy];
187  }
188  }
189 
190  }
191 
192 }
193 /**
194  * Update boundary condition for current time step
195  * \param iteration - iteration number
196  * \param PSD_2D_Slice - 2d slice of PSD for that BC
197  */
198 void BoundaryCondition::Update(int iteration, Matrix2D<double> PSD_2D_Slice, double time, double dt) {
199  int ix, iy;
200 
201  if (type == "BCT_MULTIPLE_FILES") {
202  if (time == -1 || dt == -1)
203  throw error_msg("BC_ERR", "To use BCT_MULTIPLE_FILES, time and dt need to be specified when the function is called.");
204  // Find the next time to change BC
205  int i;
206  for (i = 0; i < BC_change_time.size() && BC_change_time[i] < time; i++) {};
207  // Check if it's time to load next BC, i.e. the time to change is between the previous and current time steps
208  if (time - dt < BC_change_time[i] && time >= BC_change_time[i]) {
209  // Load next BC
210  ifstream input(BC_filename[i].c_str());
211  // checking if file was opened
212  if (input == NULL) {
213  throw error_msg("BC_FILE_OPEN_ERR", "Boundary conditions file %s not found.", BC_filename[i].c_str());
214  }
215  Output::echo(3, "Updating a boundary condition: %s", BC_filename[i].c_str());
216  // Skip header of the file
217  string inBuf;
218  input >> inBuf;
219  while (strcasecmp(inBuf.c_str(), "ZONE") != 0 && !input.eof() ) input >> inBuf;
220  input.ignore(9999, '\n'); // read to the end of the line with 'zone'
221  double tmp;
222  for (ix = 0; ix < arr.size_x; ix++) {
223  for (iy = 0; iy < arr.size_y; iy++) {
224  input >> tmp;
225  arr[ix][iy] = tmp;
226  }
227  }
228  }
229  }
230 
231 }