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(file_name, delimiter),
29 "Error in reading time data");
30}
31
33 std::string delimiter) {
35 MOFEM_LOG_CHANNEL("WORLD");
36 // Set the argument found flag to false as default
37 PetscBool arg_found = PETSC_FALSE;
38 char time_file_name[255] = {'\0'};
39 // Check to see if a filename has been provided
40 if (fileName.empty()) {
41 // If no filename, look for command line argument
42 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, fileNameFlag.c_str(),
43 time_file_name, 255, &arg_found);
44 if (arg_found) {
45 fileName = std::string(time_file_name);
46 }
47 } else {
48 // Set the command line flag to true for correct flow control using provided
49 // filename
50 arg_found = PETSC_TRUE;
51 }
52 if (!arg_found && fileName.empty() && errorIfFileNotGiven) {
53 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
54 "*** ERROR %s (time_data FILE NEEDED)", fileName.c_str());
55 } else if (!arg_found && fileName.empty() && !errorIfFileNotGiven) {
56 MOFEM_LOG_C("WORLD", Sev::warning,
57 "The %s file not provided. Loading scaled with time.",
58 fileName.c_str());
59 scalingMethod = [this](double time) { return this->getLinearScale(time); };
61 }
62
63 std::ifstream in_file_stream(fileName);
64 if (!in_file_stream.is_open() && errorIfFileNotGiven) {
65 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
66 "*** ERROR data file < %s > open unsuccessful", fileName.c_str());
67 } else if (!in_file_stream.is_open() && !errorIfFileNotGiven) {
68 MOFEM_LOG("WORLD", Sev::warning)
69 << "*** Warning data file " << fileName
70 << " open unsuccessful. Using linear time scaling.";
71 scalingMethod = [this](double time) { return this->getLinearScale(time); };
73 }
74
75 in_file_stream.seekg(0);
76 std::string line;
77 double time = 0.0, value = 0.0;
78 tSeries[time] = value;
79
80 std::regex rgx(delimiter.c_str());
81 std::sregex_token_iterator end;
82 while (std::getline(in_file_stream, line)) {
83 std::sregex_token_iterator iter(line.begin(), line.end(), rgx, -1);
84 auto time_str = iter;
85 auto value_str = ++iter;
86 if (time_str == end || value_str == end) {
87 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
88 "*** ERROR read data file error (check input time data file) ");
89 }
90 time = std::stod(time_str->str());
91 value = std::stod(value_str->str());
92 tSeries[time] = value;
93 }
94 in_file_stream.close();
95 scalingMethod = [this](double time) { return this->getScaleFromData(time); };
96 if (in_file_stream.is_open()) {
97 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
98 "*** ERROR file close unsuccessful");
99 }
101}
102
103double TimeScale::getLinearScale(const double time) { return time; }
104
105double TimeScale::getScaleFromData(const double time) {
106 if (tSeries.empty())
107 return 0.;
108
109 auto it = tSeries.find(time);
110 if (it != tSeries.end()) {
111 return it->second;
112 } else {
113 auto it = tSeries.lower_bound(time);
114 if (it == tSeries.end()) {
115 return (--it)->second;
116 }
117 if (it == tSeries.begin()) {
118 return it->second;
119 }
120 auto upper = *(it);
121 it--;
122 if (it == tSeries.end()) {
123 return upper.second;
124 }
125 auto lower = *it;
126 double t = (time - lower.first) / (upper.first - lower.first);
127 double scale1 = upper.second;
128 double scale0 = lower.second;
129 return scale0 + t * (scale1 - scale0);
130 }
131}
132
133double TimeScale::getScale(const double time) { return scalingMethod(time); }
134
135// Vector scale
136
137// template <int SPACE_DIM> TimeScaleVector<SPACE_DIM>::TimeScaleVector() {}
138template <int SPACE_DIM>
140 bool error_if_file_not_given)
141 : readFile(0), debug(0), nAme(name),
142 errorIfFileNotGiven(error_if_file_not_given) {
143 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
144}
145
146template <int SPACE_DIM>
148 bool error_if_file_not_given)
149 : readFile(0), debug(0), errorIfFileNotGiven(error_if_file_not_given) {
150 nAme = name + std::to_string(ms_id);
151 CHK_THROW_MESSAGE(timeData(), "Error in reading time data");
152}
153
156
157 char time_file_name[255];
158 PetscBool flg = PETSC_FALSE;
159 if (!nAme.empty())
160 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
161 time_file_name, 255, &flg);
162
163 if (!flg) {
164 if (errorIfFileNotGiven)
165 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
166 "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
168 }
169
170 FILE *time_data = fopen(time_file_name, "r");
171 if (time_data == NULL) {
172 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
173 "*** ERROR data file < %s > open unsuccessful", time_file_name);
174 }
175 double no1 = 0.0;
176 FTensor::Index<'i', 3> i;
177 (tSeries[no1])(i) = 0.;
178 while (!feof(time_data)) {
179 FTensor::Tensor1<double, 3> no2{0., 0., 0.};
180 int n =
181 fscanf(time_data, "%lf %lf %lf %lf", &no1, &no2(0), &no2(1), &no2(2));
182 if (n < 0) {
183 fgetc(time_data);
184 continue;
185 }
186 if (n != 4) {
187 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
188 "*** ERROR read data file error (check input time data file) "
189 "{ n = %d }",
190 n);
191 }
192 (tSeries[no1])(i) = no2(i);
193 }
194 int r = fclose(time_data);
195
196 MOFEM_LOG_CHANNEL("WORLD");
197 for (auto &[ts, vec] : tSeries) {
198 MOFEM_TAG_AND_LOG_C("WORLD", Sev::verbose, "TimeScaleVector",
199 "** read vector %3.2e time %3.2e %3.2e %3.2e", ts,
200 vec(0), vec(1), vec(2));
201 }
202 MOFEM_LOG_CHANNEL("WORLD");
203
204 if (r != 0) {
205 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
206 "*** ERROR file close unsuccessful");
207 }
208 readFile = 1;
209
210 if (readFile == 1)
211 scalingMethod = [this](double time) {
212 return this->getVectorFromData(time);
213 };
214
216}
217
218template <int SPACE_DIM>
221 return scalingMethod(time);
222}
223
224template <int SPACE_DIM>
227
230 FTensor::Tensor1<double, 3> acc0 = tSeries.begin()->second;
231
234
235 double t0 = 0, t1, dt;
236 for (auto &[ts, vec] : tSeries) {
237 if (ts > time) {
238 t1 = ts;
239 dt = time - t0;
240 acc(I) = acc0(I) + ((vec(I) - acc0(I)) / (t1 - t0)) * dt;
241 break;
242 }
243 t0 = ts;
244 acc0 = vec;
245 acc = acc0;
246 }
247 Nf(i) = acc(i);
248 return Nf;
249}
250
251template class TimeScaleVector<3>;
252template class TimeScaleVector<2>;
253
254} // 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.
static const bool debug
FTensor::Index< 'n', SPACE_DIM > n
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
double dt
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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 fileName, 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)