v0.13.0
TimeForceScale.hpp
Go to the documentation of this file.
1 /* \brief TimeForceScale.hpp
2  *
3  * Edited and modified by Hassan.
4  *
5  * This is not exactly procedure for linear elastic dynamics, since jacobian is
6  * evaluated at every time step and SNES procedure is involved. However it is
7  * implemented like that, to test methodology for general nonlinear problem.
8  *
9  */
10 
11 /* This file is part of MoFEM.
12  * MoFEM is free software: you can redistribute it and/or modify it under
13  * the terms of the GNU Lesser General Public License as published by the
14  * Free Software Foundation, either version 3 of the License, or (at your
15  * option) any later version.
16  *
17  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
18  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
19  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
20  * License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public
23  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
24 
25 #ifndef __TIMEFORCESCALE_HPP__
26 #define __TIMEFORCESCALE_HPP__
27 
28 /** \brief Force scale operator for reading two columns
29  */
31 
32  std::map<double, double> tSeries;
34  string nAme;
36 
37  TimeForceScale(string name = "-my_time_data_file",
38  bool error_if_file_not_given = false)
39  : readFile(0), debug(0), nAme(name),
40  errorIfFileNotGiven(error_if_file_not_given) {
41 
42  ierr = timeData();
43  CHKERRABORT(PETSC_COMM_WORLD, ierr);
44  }
45 
46  PetscBool fLg;
47 
50  char time_file_name[255];
51  ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
52  time_file_name, 255, &fLg);
53  CHKERRG(ierr);
54  if (!fLg && errorIfFileNotGiven) {
55  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
56  "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
57  }
58  if (!fLg) {
59  MOFEM_LOG_C("WORLD", Sev::warning,
60  "The %s file not provided. Loading scaled with time.",
61  nAme.c_str());
63  }
64  FILE *time_data = fopen(time_file_name, "r");
65  if (time_data == NULL) {
66  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
67  "*** ERROR data file < %s > open unsuccessful", time_file_name);
68  }
69  double no1 = 0.0, no2 = 0.0;
70  tSeries[no1] = no2;
71  while (!feof(time_data)) {
72  int n = fscanf(time_data, "%lf %lf", &no1, &no2);
73  if ((n <= 0) || ((no1 == 0) && (no2 == 0))) {
74  fgetc(time_data);
75  continue;
76  }
77  if (n != 2) {
78  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79  "*** ERROR read data file error (check input time data file) "
80  "{ n = %d }",
81  n);
82  }
83  tSeries[no1] = no2;
84  }
85  int r = fclose(time_data);
86  if (debug) {
87  std::map<double, double>::iterator tit = tSeries.begin();
88  for (; tit != tSeries.end(); tit++) {
89  PetscPrintf(PETSC_COMM_WORLD, "** read time series %3.2e time %3.2e\n",
90  tit->first, tit->second);
91  }
92  }
93  if (r != 0) {
94  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
95  "*** ERROR file close unsuccessful");
96  }
97  readFile = 1;
99  }
100 
101  MoFEMErrorCode getForceScale(const double ts_t, double &scale) {
103  if (!fLg) {
104  scale = ts_t; // scale with time, by default
106  }
107  if (readFile == 0) {
108  SETERRQ(PETSC_COMM_SELF, 1, "data file not read");
109  }
110  scale = 0;
111  double t0 = 0, t1, s0 = tSeries[0], s1, dt;
112  std::map<double, double>::iterator tit = tSeries.begin();
113  for (; tit != tSeries.end(); tit++) {
114  if (tit->first > ts_t) {
115  t1 = tit->first;
116  s1 = tit->second;
117  dt = ts_t - t0;
118  scale = s0 + ((s1 - s0) / (t1 - t0)) * dt;
119  break;
120  }
121  t0 = tit->first;
122  s0 = tit->second;
123  scale = s0;
124  }
126  }
127 
128  /**
129  * @brief Scale force the right hand vector
130  *
131  * @param fe
132  * @param Nf
133  * @return MoFEMErrorCode
134  */
137  double scale;
138  const double ts_t = fe->ts_t;
139  CHKERR getForceScale(ts_t, scale);
140  Nf *= scale;
142  }
143 };
144 
146 
147  std::map<double, VectorDouble> tSeries;
149  string nAme;
150 
151  TimeAccelerogram(string name = "-my_accelerogram")
152  : readFile(0), debug(0), nAme(name) {
153 
154  ierr = timeData();
155  CHKERRABORT(PETSC_COMM_WORLD, ierr);
156  }
157 
160  char time_file_name[255];
161  PetscBool flg = PETSC_TRUE;
162  ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
163  time_file_name, 255, &flg);
164  CHKERRG(ierr);
165  if (flg != PETSC_TRUE) {
166  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
167  "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
168  }
169  FILE *time_data = fopen(time_file_name, "r");
170  if (time_data == NULL) {
171  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
172  "*** ERROR data file < %s > open unsuccessful", time_file_name);
173  }
174  double no1 = 0.0;
175  VectorDouble no2(3);
176  tSeries[no1] = no2;
177  while (!feof(time_data)) {
178  int n =
179  fscanf(time_data, "%lf %lf %lf %lf", &no1, &no2[0], &no2[1], &no2[2]);
180  if (n < 0) {
181  fgetc(time_data);
182  continue;
183  }
184  if (n != 4) {
185  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
186  "*** ERROR read data file error (check input time data file) "
187  "{ n = %d }",
188  n);
189  }
190  tSeries[no1] = no2;
191  }
192  int r = fclose(time_data);
193  if (debug) {
194  std::map<double, VectorDouble>::iterator tit = tSeries.begin();
195  for (; tit != tSeries.end(); tit++) {
196  PetscPrintf(PETSC_COMM_WORLD,
197  "** read accelerogram %3.2e time %3.2e %3.2e %3.2e\n",
198  tit->first, tit->second[0], tit->second[1], tit->second[2]);
199  }
200  }
201  if (r != 0) {
202  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR file close unsuccessful");
203  }
204  readFile = 1;
206  }
207 
210  if (readFile == 0) {
211  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data file not read");
212  }
213  double ts_t = fe->ts_t;
214  VectorDouble acc(3);
215  VectorDouble acc0 = tSeries[0], acc1(3);
216  double t0 = 0, t1, dt;
217  std::map<double, VectorDouble>::iterator tit = tSeries.begin();
218  for (; tit != tSeries.end(); tit++) {
219  if (tit->first > ts_t) {
220  t1 = tit->first;
221  acc1 = tit->second;
222  dt = ts_t - t0;
223  acc = acc0 + ((acc1 - acc0) / (t1 - t0)) * dt;
224  break;
225  }
226  t0 = tit->first;
227  acc0 = tit->second;
228  acc = acc0;
229  }
230  Nf += acc;
232  }
233 };
234 
235 #endif // __TIMEFORCESCALE_HPP__
static Index< 'n', 3 > n
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
double dt
Definition: heat_method.cpp:40
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
const double r
rate factor
double scale
Definition: plastic.cpp:103
Class used to scale loads, f.e. in arc-length control.
std::map< double, VectorDouble > tSeries
TimeAccelerogram(string name="-my_accelerogram")
MoFEMErrorCode timeData()
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &Nf)
Force scale operator for reading two columns.
TimeForceScale(string name="-my_time_data_file", bool error_if_file_not_given=false)
std::map< double, double > tSeries
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &Nf)
Scale force the right hand vector.
MoFEMErrorCode timeData()
MoFEMErrorCode getForceScale(const double ts_t, double &scale)