00001
00008 #include "BoundaryConditions.h"
00009 #include <math.h>
00010
00011 BoundaryCondition::BoundaryCondition(int time_size, int size_two, int size_three, ParamStructure::BoundaryCondition parameters) : Matrix3D<double>(time_size, size_two, size_three) {
00012
00013 Initialize(parameters);
00014 }
00015
00024 void BoundaryCondition::Initialize(ParamStructure::BoundaryCondition parameters) {
00025
00026
00027 if (!Matrix3D<double>::initialized) {
00028 throw error_msg("MEMORY_IS_NOT_ALLOCATED", "Boundary conditions memory is not allocated");
00029 }
00030
00031 BC_parameters = parameters;
00032 this->initialType = parameters.type;
00033
00034 if (parameters.type == "BCT_CONSTANT_VALUE") {
00035
00036 this->calculationType = "BCT_CONSTANT_VALUE";
00037 this->value = parameters.value;
00038
00039
00040
00041 (*this) = parameters.value;
00042 } else if (parameters.type == "BCT_CONSTANT_DERIVATIVE") {
00043
00044 this->calculationType = "BCT_CONSTANT_DERIVATIVE";
00045 this->value = parameters.value;
00046
00047
00048
00049 (*this) = parameters.value;
00050 } else if (parameters.type == "BCT_CONSTANT_PSD") {
00051
00052 this->calculationType = "BCT_CONSTANT_VALUE";
00053 this->value = parameters.value;
00054 } else if (parameters.type == "BCT_FILE" || parameters.type == "BCT_FILE_GRID") {
00055
00056 this->calculationType = "BCT_CONSTANT_VALUE";
00057 this->filename = parameters.filename;
00058 } else if (parameters.type == "BCT_FILE_CHANGES") {
00059
00060 this->calculationType = "BCT_CONSTANT_VALUE";
00061 this->filename = parameters.filename;
00062 } else {
00063 throw error_msg("BOUNDARY_CONDITION_TYPE_ERROR", "Unknown boundary condition type");
00064 }
00065
00066 initialized = true;
00067 }
00068
00078 void BoundaryCondition::MakeBoundaryCondition(Matrix2D<double> psd2DSlice, Matrix2D<double> gridElement1, Matrix2D<double> gridElement2) {
00079 int it, iy, iz;
00080
00081 if (!initialized) {
00082 throw error_msg("BOUNDARY_CONDITION_INITIALIZATION_ERROR", "Boundary condition was not initialized");
00083 }
00084 if (initialType == "BCT_CONSTANT_VALUE") {
00085
00086
00087 (*this) = VF::max(value, VC::zero_f);
00088
00089 this->calculationType = "BCT_CONSTANT_VALUE";
00090 } else if (initialType == "BCT_CONSTANT_PSD") {
00091
00092 for (it = 0; it < this->size_x; it++) {
00093 for (iy = 0; iy < this->size_y; iy++) {
00094 for (iz = 0; iz < this->size_z; iz++) {
00095 (*this)[it][iy][iz] = VF::max(psd2DSlice[iy][iz], VC::zero_f);
00096 }
00097 }
00098 }
00099 this->calculationType = "BCT_CONSTANT_VALUE";
00100 } else if (initialType == "BCT_CONSTANT_DERIVATIVE") {
00101 (*this) = value;
00102 this->calculationType = "BCT_CONSTANT_DERIVATIVE";
00103 } else if (initialType == "BCT_FILE" || initialType == "BCT_FILE_GRID" || initialType == "BCT_FILE_CHANGES") {
00104
00105 } else {
00106 throw error_msg("BCT_ERROR", "Unknown boundary condition type: %d", initialType.c_str());
00107 }
00108 }
00109
00118 void BoundaryCondition::LoadBoundaryCondition(Matrix2D<double> gridElement1, Matrix2D<double> gridElement2) {
00119 int it, iy, iz;
00120 double err = 1.e-3;
00121
00122 if (!initialized) {
00123 throw error_msg("BOUNDARY_CONDITION_INITIALIZATION_ERROR", "Boundary condition was not initialized");
00124 }
00125 if (initialType == "BCT_FILE" || initialType == "BCT_FILE_GRID" || initialType == "BCT_FILE_CHANGES") {
00126 Output::echo(0, "Loading boundary conditions from files... ");
00127
00128
00129 {
00130 ifstream input(filename.c_str());
00131
00132 if (input == NULL) {
00133 throw error_msg("BC_FILE_OPEN_ERR", "Boundary conditions file %s not found.", filename.c_str());
00134 }
00135
00136 input.ignore(9999, '\n');
00137 input.ignore(9999, '\n');
00138 double loadedTime, loadedAx1, loadedAx2, loadedValue;
00139
00140 for (it = 0; it < this->size_x; it++) {
00141 for (iy = 0; iy < this->size_y; iy++) {
00142 for (iz = 0; iz < this->size_z; iz++) {
00143 if (initialType == "BCT_FILE_GRID") input >> loadedTime >> loadedAx1 >> loadedAx2;
00144 input >> loadedValue;
00145
00146 if (initialType == "BCT_FILE_GRID" && (fabs(loadedAx1 - gridElement1[iy][iz]) > err || fabs(loadedAx2 - gridElement2[iy][iz]) > err)) {
00147 if (input.eof()) {
00148 throw error_msg("BC_LOAD_GRID_ERR", "Loading %s: not enough time steps in the boundary condition file (%d found/%d needed).\n", filename.c_str(), it, this->size_x);
00149 } else {
00150 throw error_msg("BC_LOAD_GRID_ERR", "Loading %s: grid mismatch.\nLoaded: ax1=%e, ax2=%e\nGrid: ax1=%e, ax2=%e\n(it=%d, i1=%d, i2=%d)\n", filename.c_str(), loadedAx1, loadedAx2, gridElement1[iy][iz], gridElement2[iy][iz], it, iy, iz);
00151 }
00152 } else {
00153 (*this)[it][iy][iz] = loadedValue;
00154 }
00155
00156 input.ignore(9999, '\n');
00157 }
00158 }
00159 }
00160 input.close();
00161 }
00162 (*this).writeToFile("./Debug/" + (*this).name + "_BC.dat");
00163 Output::echo(0, "done.\n");
00164 }
00165 }
00171 void BoundaryCondition::Update(int iteration, Matrix2D<double> PSD_2D_Slice) {
00172 int iy, iz;
00173 if (initialType == "BCT_FILE_CHANGES") {
00174 Output::echo(0, "Changing boundary conditions %s.\n", PSD_2D_Slice.name.c_str());
00175
00176 for (iy = 0; iy < this->size_y; iy++) {
00177 for (iz = 0; iz < this->size_z; iz++) {
00178
00179
00180
00181 (*this)[iteration][iy][iz] = PSD_2D_Slice[iy][iz];
00182 }
00183 }
00184 }
00185 }
00186
00187
00193 BoundaryCondition& BoundaryCondition::operator= (double val) {
00194 Matrix3D<double>::operator=(val);
00195 return *this;
00196 }
00197
00203
00204 BoundaryCondition& BoundaryCondition::operator= (Matrix3D<double> M) {
00205 Matrix3D<double>::operator=(M);
00206 return *this;
00207 }
00208
00214 BoundaryCondition BoundaryCondition::operator* (double val) {
00215 BoundaryCondition tmp(*this);
00216 tmp = tmp.Matrix3D<double>::operator* (val);
00217 return tmp;
00218 }