v0.14.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 : 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 MOFEM_LOG_CHANNEL("WORLD");
34 PetscBool arg_found = PETSC_FALSE;
35 char time_file_name[255] = {'\0'};
36 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, fileNameFlag.c_str(),
37 time_file_name, 255, &arg_found);
38 if (arg_found) {
39 fileName = std::string(time_file_name);
40 }
41 if (!arg_found && fileName.empty() && errorIfFileNotGiven) {
42 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
43 "*** ERROR %s (time_data FILE NEEDED)", fileName.c_str());
44 } else if (!arg_found && fileName.empty() && !errorIfFileNotGiven) {
45 MOFEM_LOG_C("WORLD", Sev::warning,
46 "The %s file not provided. Loading scaled with time.",
47 fileName.c_str());
48 scalingMethod = [this](double time) { return this->getLinearScale(time); };
50 }
51
52 std::ifstream in_file_stream(fileName);
53 if (!in_file_stream.is_open() && errorIfFileNotGiven) {
54 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
55 "*** ERROR data file < %s > open unsuccessful", fileName.c_str());
56 } else if (!in_file_stream.is_open() && !errorIfFileNotGiven) {
57 MOFEM_LOG("WORLD", Sev::warning)
58 << "*** Warning data file " << fileName
59 << " open unsuccessful. Using linear time scaling.";
60 scalingMethod = [this](double time) { return this->getLinearScale(time); };
62 }
63
64 in_file_stream.seekg(0);
65 std::string line;
66 double time = 0.0, value = 0.0;
67 tSeries[time] = value;
68
69 std::regex rgx(delimiter.c_str());
70 std::sregex_token_iterator end;
71 while (std::getline(in_file_stream, line)) {
72 std::sregex_token_iterator iter(line.begin(), line.end(), rgx, -1);
73 auto time_str = iter;
74 auto value_str = ++iter;
75 if (time_str == end || value_str == end) {
76 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
77 "*** ERROR read data file error (check input time data file) ");
78 }
79 time = std::stod(time_str->str());
80 value = std::stod(value_str->str());
81 tSeries[time] = value;
82 }
83 in_file_stream.close();
84 scalingMethod = [this](double time) { return this->getScaleFromData(time); };
85 if (in_file_stream.is_open()) {
86 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
87 "*** ERROR file close unsuccessful");
88 }
90}
91
92double TimeScale::getLinearScale(const double time) { return time; }
93
94double TimeScale::getScaleFromData(const double time) {
95 if (tSeries.empty())
96 return 0.;
97
98 auto it = tSeries.find(time);
99 if (it != tSeries.end()) {
100 return it->second;
101 } else {
102 auto it = tSeries.lower_bound(time);
103 if (it == tSeries.end()) {
104 return (--it)->second;
105 }
106 if (it == tSeries.begin()) {
107 return it->second;
108 }
109 auto upper = *(it);
110 it--;
111 if (it == tSeries.end()) {
112 return upper.second;
113 }
114 auto lower = *it;
115 double t = (time - lower.first) / (upper.first - lower.first);
116 double scale1 = upper.second;
117 double scale0 = lower.second;
118 return scale0 + t * (scale1 - scale0);
119 }
120}
121
122double TimeScale::getScale(const double time) { return scalingMethod(time); }
123
124// Vector scale
125
126// template <int SPACE_DIM> TimeScaleVector<SPACE_DIM>::TimeScaleVector() {}
127template <int SPACE_DIM>
129 bool error_if_file_not_given)
130 : readFile(0), debug(0), nAme(name),
131 errorIfFileNotGiven(error_if_file_not_given) {
132 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
133}
134
135template <int SPACE_DIM>
137 bool error_if_file_not_given)
138 : readFile(0), debug(0), errorIfFileNotGiven(error_if_file_not_given) {
139 nAme = name + std::to_string(ms_id);
140 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
141}
142
145
146 char time_file_name[255];
147 PetscBool flg = PETSC_FALSE;
148 if (!nAme.empty())
149 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
150 time_file_name, 255, &flg);
151
152 if (!flg) {
153 if (errorIfFileNotGiven)
154 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
155 "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
157 }
158
159 FILE *time_data = fopen(time_file_name, "r");
160 if (time_data == NULL) {
161 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
162 "*** ERROR data file < %s > open unsuccessful", time_file_name);
163 }
164 double no1 = 0.0;
166 (tSeries[no1])(i) = 0.;
167 while (!feof(time_data)) {
168 FTensor::Tensor1<double, 3> no2{0., 0., 0.};
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])(i) = no2(i);
182 }
183 int r = fclose(time_data);
184
185 MOFEM_LOG_CHANNEL("WORLD");
186 for (auto &[ts, vec] : tSeries) {
187 MOFEM_TAG_AND_LOG_C("WORLD", Sev::verbose, "TimeScaleVector",
188 "** read vector %3.2e time %3.2e %3.2e %3.2e", ts,
189 vec(0), vec(1), vec(2));
190 }
191 MOFEM_LOG_CHANNEL("WORLD");
192
193 if (r != 0) {
194 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
195 "*** ERROR file close unsuccessful");
196 }
197 readFile = 1;
198
199 if (readFile == 1)
200 scalingMethod = [this](double time) {
201 return this->getVectorFromData(time);
202 };
203
205}
206
207template <int SPACE_DIM>
210 return scalingMethod(time);
211}
212
213template <int SPACE_DIM>
216
219 FTensor::Tensor1<double, 3> acc0 = tSeries.begin()->second;
220
223
224 double t0 = 0, t1, dt;
225 for (auto &[ts, vec] : tSeries) {
226 if (ts > time) {
227 t1 = ts;
228 dt = time - t0;
229 acc(I) = acc0(I) + ((vec(I) - acc0(I)) / (t1 - t0)) * dt;
230 break;
231 }
232 t0 = ts;
233 acc0 = vec;
234 acc = acc0;
235 }
236 Nf(i) = acc(i);
237 return Nf;
238}
239
240template class TimeScaleVector<3>;
241template class TimeScaleVector<2>;
242
243} // namespace MoFEM
#define MOFEM_TAG_AND_LOG_C(channel, severity, tag, format,...)
Tag and log in channel.
Definition: LogManager.hpp:370
#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
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#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 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
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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 IntegrationType I
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)