v0.13.2
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 : TimeScale(file_name, defaultDelimiter, error_if_file_not_given) {}
24
25TimeScale::TimeScale(std::string file_name, std::string delimiter,
26 bool error_if_file_not_given)
27 : fileName(file_name), errorIfFileNotGiven(error_if_file_not_given) {
28 CHK_THROW_MESSAGE(timeData(delimiter), "Error in reading time data");
29}
30
31MoFEMErrorCode TimeScale::timeData(std::string delimiter) {
33 PetscBool arg_found = PETSC_FALSE;
34 char time_file_name[255] = {'\0'};
35 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, fileNameFlag.c_str(),
36 time_file_name, 255, &arg_found);
37 if (arg_found) {
38 fileName = std::string(time_file_name);
39 }
40 if (!arg_found && fileName.empty() && errorIfFileNotGiven) {
41 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
42 "*** ERROR %s (time_data FILE NEEDED)", fileName.c_str());
43 } else if (!arg_found && fileName.empty() && !errorIfFileNotGiven) {
44 MOFEM_LOG_C("WORLD", Sev::warning,
45 "The %s file not provided. Loading scaled with time.",
46 fileName.c_str());
47 scalingMethod = [this](double time) { return this->getLinearScale(time); };
49 }
50
51 std::ifstream in_file_stream(fileName);
52 if (!in_file_stream.is_open() && errorIfFileNotGiven) {
53 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
54 "*** ERROR data file < %s > open unsuccessful", fileName.c_str());
55 } else if (!in_file_stream.is_open() && !errorIfFileNotGiven) {
56 MOFEM_LOG("WORLD", Sev::warning)
57 << "*** Warning data file " << fileName
58 << " open unsuccessful. Using linear time scaling.";
59 scalingMethod = [this](double time) { return this->getLinearScale(time); };
60 }
61 in_file_stream.seekg(0);
62 std::string line;
63 double time = 0.0, value = 0.0;
64 tSeries[time] = value;
65
66 std::regex rgx(delimiter.c_str());
67 std::sregex_token_iterator end;
68 while (std::getline(in_file_stream, line)) {
69 std::sregex_token_iterator iter(line.begin(), line.end(), rgx, -1);
70 auto time_str = iter;
71 auto value_str = ++iter;
72 if (time_str == end || value_str == end) {
73 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
74 "*** ERROR read data file error (check input time data file) ");
75 }
76 time = std::stod(time_str->str());
77 value = std::stod(value_str->str());
78 tSeries[time] = value;
79 }
80 in_file_stream.close();
81 scalingMethod = [this](double time) { return this->getScaleFromData(time); };
82 if (in_file_stream.is_open()) {
83 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
84 "*** ERROR file close unsuccessful");
85 }
87}
88
89double TimeScale::getLinearScale(const double time) { return time; }
90
91double TimeScale::getScaleFromData(const double time) {
92 if (tSeries.empty())
93 return 0.;
94
95 auto it = tSeries.find(time);
96 if (it != tSeries.end()) {
97 return it->second;
98 } else {
99 auto it = tSeries.lower_bound(time);
100 if (it == tSeries.end()) {
101 return (--it)->second;
102 }
103 if (it == tSeries.begin()) {
104 return it->second;
105 }
106 auto upper = *(it);
107 it--;
108 if (it == tSeries.end()) {
109 return upper.second;
110 }
111 auto lower = *it;
112 double t = (time - lower.first) / (upper.first - lower.first);
113 double scale1 = upper.second;
114 double scale0 = lower.second;
115 return scale0 + t * (scale1 - scale0);
116 }
117}
118
119double TimeScale::getScale(const double time) { return scalingMethod(time); }
120
121// Vector scale
122
123// template <int SPACE_DIM> TimeScaleVector<SPACE_DIM>::TimeScaleVector() {}
124template <int SPACE_DIM>
126 bool error_if_file_not_given)
127 : readFile(0), debug(0), nAme(name),
128 errorIfFileNotGiven(error_if_file_not_given) {
129 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
130}
131
132template <int SPACE_DIM>
134 bool error_if_file_not_given) : readFile(0), debug(0),
135 errorIfFileNotGiven(error_if_file_not_given) {
136 nAme = name + std::to_string(ms_id);
137 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
138}
139
140template <int SPACE_DIM>
142
144
145 char time_file_name[255];
146 PetscBool flg = PETSC_FALSE;
147 if (!nAme.empty())
148 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
149 time_file_name, 255, &flg);
150
151 if (!flg) {
152 if (errorIfFileNotGiven)
153 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
154 "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
156 }
157
158 FILE *time_data = fopen(time_file_name, "r");
159 if (time_data == NULL) {
160 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
161 "*** ERROR data file < %s > open unsuccessful", time_file_name);
162 }
163 double no1 = 0.0;
166 no2(i) = 0.;
167 tSeries[no1] = no2;
168 while (!feof(time_data)) {
169 int n =
170 fscanf(time_data, "%lf %lf %lf %lf", &no1, &no2(0), &no2(1), &no2(2));
171 if (n < 0) {
172 fgetc(time_data);
173 continue;
174 }
175 if (n != 4) {
176 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
177 "*** ERROR read data file error (check input time data file) "
178 "{ n = %d }",
179 n);
180 }
181 tSeries[no1] = no2;
182 }
183 int r = fclose(time_data);
184 if (debug) {
185
186 for (auto &[ts, vec] : tSeries) {
187 PetscPrintf(PETSC_COMM_WORLD,
188 "** read vector %3.2e time %3.2e %3.2e %3.2e\n",
189 ts, vec(0), vec(1), vec(2));
190 }
191 }
192 if (r != 0) {
193 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR file close unsuccessful");
194 }
195 readFile = 1;
196
197 if (readFile == 1)
198 scalingMethod = [this](double time) {
199 return this->getVectorFromData(time);
200 };
201
203}
204
205
206template <int SPACE_DIM>
209 return scalingMethod(time);
210}
211
212template <int SPACE_DIM>
215
216 // if (readFile == 0) {
217 // CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Data file not read");
218 // }
219
221 FTensor::Tensor1<double, SPACE_DIM> acc0 = tSeries.begin()->second;
224
225 double t0 = 0, t1, dt;
226 for (auto &[ts, vec] : tSeries) {
227 if (ts > time) {
228 t1 = ts;
229 dt = time - t0;
230 acc(i) = acc0(i) + ((vec(i) - acc0(i)) / (t1 - t0)) * dt;
231 break;
232 }
233 t0 = ts;
234 acc0 = vec;
235 acc = acc0;
236 }
237 Nf(i) = acc(i);
238 return Nf;
239}
240
241template class TimeScaleVector<3>;
242template class TimeScaleVector<2>;
243
244} // namespace MoFEM
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
static const bool debug
FTensor::Index< 'n', SPACE_DIM > n
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FTensor::Index< 'i', SPACE_DIM > i
double dt
Definition: heat_method.cpp:26
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
constexpr double t
plate stiffness
Definition: plate.cpp:59
virtual double getScale(const double time)
Get scaling at given time.
Force scale operator for reading two columns.
std::string fileName
double getLinearScale(const double time)
Returns the value of time.
double getScale(const double time)
Get scaling at a given time.
double getScaleFromData(const double time)
Get scaling at a given time when the scalar values have been provided. Uses linear interpolation on t...
TimeScale(std::string file_name="", bool error_if_file_not_given=false)
TimeScale constructor.
std::string fileNameFlag
static const std::string defaultDelimiter
comma or space
std::function< double(double)> scalingMethod
MoFEMErrorCode timeData(std::string delimiter)
std::map< double, double > tSeries
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)