v0.15.0
Loading...
Searching...
No Matches
ScalingMethod.cpp
Go to the documentation of this file.
1/**
2 * @file ScalingMethod.cpp
3 * @brief
4 * @version 0.1
5 * @date 2022-08-12
6 *
7 * @copyright Copyright (c) 2022
8 *
9 */
10
11#include <MoFEM.hpp>
12
13namespace MoFEM {
14
15double ScalingMethod::getScale(const double time) {
16 THROW_MESSAGE("getScale not implemented");
17}
18
19const std::string TimeScale::defaultDelimiter =
20 "(\\s*,\\s*|\\s+)"; ///< comma or space
21
22TimeScale::TimeScale(std::string file_name, bool error_if_file_not_given,
23 ScalingFun def_scaling_fun)
24 : TimeScale(file_name, defaultDelimiter, error_if_file_not_given,
25 def_scaling_fun) {}
26
27TimeScale::TimeScale(std::string file_name, std::string delimiter,
28 bool error_if_file_not_given, ScalingFun def_scaling_fun)
29 : fileName(file_name), errorIfFileNotGiven(error_if_file_not_given),
30 defScalingMethod(def_scaling_fun) {
31 CHK_THROW_MESSAGE(timeData(file_name, delimiter),
32 "Error in reading time data");
33}
34
36 std::string delimiter) {
38 MOFEM_LOG_CHANNEL("WORLD");
39 // Set the argument found flag to false as default
40 char time_file_name[255] = {'\0'};
41 // Check to see if a filename has been provided
42 if (fileName.empty()) {
43 // If no filename, look for command line argument
44 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, fileNameFlag.c_str(),
45 time_file_name, 255, &argFileScale);
46 if (argFileScale) {
47 fileName = std::string(time_file_name);
48 }
49 } else {
50 // Set the command line flag to true for correct flow control using provided
51 // filename
52 argFileScale = PETSC_TRUE;
53 }
54 if (!argFileScale && fileName.empty() && errorIfFileNotGiven) {
55 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
56 "*** ERROR %s (time_data FILE NEEDED)", fileName.c_str());
57 } else if (!argFileScale && fileName.empty() && !errorIfFileNotGiven) {
58 MOFEM_LOG_C("WORLD", Sev::warning,
59 "The %s file not provided. Loading scaled with time.",
60 fileName.c_str());
61 scalingMethod = [this](double time) {
62 return this->defScalingMethod(time);
63 };
65 }
66
67 std::ifstream in_file_stream(fileName);
68 if (!in_file_stream.is_open() && errorIfFileNotGiven) {
69 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
70 "*** ERROR data file < %s > open unsuccessful", fileName.c_str());
71 } else if (!in_file_stream.is_open() && !errorIfFileNotGiven) {
72 MOFEM_LOG("WORLD", Sev::warning)
73 << "*** Warning data file " << fileName
74 << " open unsuccessful. Using default time scaling.";
75 scalingMethod = [this](double time) {
76 return this->defScalingMethod(time);
77 };
79 }
80
81 in_file_stream.seekg(0);
82 std::string line;
83 double time = 0.0, value = 0.0;
84 tSeries[time] = value;
85
86 std::regex rgx(delimiter.c_str());
87 std::sregex_token_iterator end;
88 while (std::getline(in_file_stream, line)) {
89 std::sregex_token_iterator iter(line.begin(), line.end(), rgx, -1);
90 auto time_str = iter;
91 auto value_str = ++iter;
92 if (time_str == end || value_str == end) {
93 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
94 "*** ERROR read data file error (check input time data file) ");
95 }
96 time = std::stod(time_str->str());
97 value = std::stod(value_str->str());
98 tSeries[time] = value;
99 }
100 in_file_stream.close();
101 scalingMethod = [this](double time) { return this->getScaleFromData(time); };
102 if (in_file_stream.is_open()) {
103 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
104 "*** ERROR file close unsuccessful");
105 }
107}
108
109double TimeScale::getScaleFromData(const double time) {
110 if (tSeries.empty())
111 return 0.;
112
113 auto it = tSeries.find(time);
114 if (it != tSeries.end()) {
115 return it->second;
116 } else {
117 auto it = tSeries.lower_bound(time);
118 if (it == tSeries.end()) {
119 return (--it)->second;
120 }
121 if (it == tSeries.begin()) {
122 return it->second;
123 }
124 auto upper = *(it);
125 it--;
126 if (it == tSeries.end()) {
127 return upper.second;
128 }
129 auto lower = *it;
130 double t = (time - lower.first) / (upper.first - lower.first);
131 double scale1 = upper.second;
132 double scale0 = lower.second;
133 return scale0 + t * (scale1 - scale0);
134 }
135}
136
137double TimeScale::getScale(const double time) { return scalingMethod(time); }
138
139// Vector scale
140
141// template <int SPACE_DIM> TimeScaleVector<SPACE_DIM>::TimeScaleVector() {}
142template <int SPACE_DIM>
144 bool error_if_file_not_given)
145 : readFile(0), debug(0), nAme(name),
146 errorIfFileNotGiven(error_if_file_not_given) {
147 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
148}
149
150template <int SPACE_DIM>
152 bool error_if_file_not_given)
153 : readFile(0), debug(0), errorIfFileNotGiven(error_if_file_not_given) {
154 nAme = name + std::to_string(ms_id);
155 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
156}
157
160
161 char time_file_name[255];
162 PetscBool flg = PETSC_FALSE;
163 if (!nAme.empty())
164 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, nAme.c_str(),
165 time_file_name, 255, &flg);
166
167 if (!flg) {
168 if (errorIfFileNotGiven)
169 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
170 "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
172 }
173
174 FILE *time_data = fopen(time_file_name, "r");
175 if (time_data == NULL) {
176 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
177 "*** ERROR data file < %s > open unsuccessful", time_file_name);
178 }
179 double no1 = 0.0;
180 FTensor::Index<'i', 3> i;
181 (tSeries[no1])(i) = 0.;
182 while (!feof(time_data)) {
183 FTensor::Tensor1<double, 3> no2{0., 0., 0.};
184 int n =
185 fscanf(time_data, "%lf %lf %lf %lf", &no1, &no2(0), &no2(1), &no2(2));
186 if (n < 0) {
187 fgetc(time_data);
188 continue;
189 }
190 if (n != 4) {
191 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
192 "*** ERROR read data file error (check input time data file) "
193 "{ n = %d }",
194 n);
195 }
196 (tSeries[no1])(i) = no2(i);
197 }
198 int r = fclose(time_data);
199
200 MOFEM_LOG_CHANNEL("WORLD");
201 for (auto &[ts, vec] : tSeries) {
202 MOFEM_TAG_AND_LOG_C("WORLD", Sev::verbose, "TimeScaleVector",
203 "** read vector %3.2e time %3.2e %3.2e %3.2e", ts,
204 vec(0), vec(1), vec(2));
205 }
206 MOFEM_LOG_CHANNEL("WORLD");
207
208 if (r != 0) {
209 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
210 "*** ERROR file close unsuccessful");
211 }
212 readFile = 1;
213
214 if (readFile == 1)
215 scalingMethod = [this](double time) {
216 return this->getVectorFromData(time);
217 };
218
220}
221
222template <int SPACE_DIM>
225 return scalingMethod(time);
226}
227
228template <int SPACE_DIM>
231
234 FTensor::Tensor1<double, 3> acc0 = tSeries.begin()->second;
235
238
239 double t0 = 0, t1, dt;
240 for (auto &[ts, vec] : tSeries) {
241 if (ts > time) {
242 t1 = ts;
243 dt = time - t0;
244 acc(I) = acc0(I) + ((vec(I) - acc0(I)) / (t1 - t0)) * dt;
245 break;
246 }
247 t0 = ts;
248 acc0 = vec;
249 acc = acc0;
250 }
251 Nf(i) = acc(i);
252 return Nf;
253}
254
255template class TimeScaleVector<3>;
256template class TimeScaleVector<2>;
257
258
259
260} // namespace MoFEM
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
#define MOFEM_LOG_C(channel, severity, format,...)
constexpr int SPACE_DIM
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
double dt
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static const bool debug
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
constexpr IntegrationType I
constexpr double t
plate stiffness
Definition plate.cpp:58
virtual double getScale(const double time)
Get scaling at given time.
Force scale operator for reading four columns (time and vector)
TimeScaleVector(std::string name="-time_vector_file", bool error_if_file_not_given=false)
MoFEMErrorCode timeData()
virtual FTensor::Tensor1< double, SPACE_DIM > getVector(const double time)
virtual FTensor::Tensor1< double, SPACE_DIM > getVectorFromData(const double time)
Force scale operator for reading two columns.
std::string fileName
ScalingFun defScalingMethod
std::function< double(double)> ScalingFun
ScalingFun scalingMethod
double getScaleFromData(const double time)
Get scaling at a given time when the scalar values have been provided. Uses linear interpolation on t...
MoFEMErrorCode timeData(std::string fileName, std::string delimiter)
std::string fileNameFlag
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
static const std::string defaultDelimiter
comma or space
std::map< double, double > tSeries