v0.9.1
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 = true)
39  : readFile(0), debug(1), 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) {
60  }
61  FILE *time_data = fopen(time_file_name, "r");
62  if (time_data == NULL) {
63  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
64  "*** ERROR data file < %s > open unsuccessful", time_file_name);
65  }
66  double no1 = 0.0, no2 = 0.0;
67  tSeries[no1] = no2;
68  while (!feof(time_data)) {
69  int n = fscanf(time_data, "%lf %lf", &no1, &no2);
70  if ((n <= 0) || ((no1 == 0) && (no2 == 0))) {
71  fgetc(time_data);
72  continue;
73  }
74  if (n != 2) {
75  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
76  "*** ERROR read data file error (check input time data file) "
77  "{ n = %d }",
78  n);
79  }
80  tSeries[no1] = no2;
81  }
82  int r = fclose(time_data);
83  if (debug) {
84  std::map<double, double>::iterator tit = tSeries.begin();
85  for (; tit != tSeries.end(); tit++) {
86  PetscPrintf(PETSC_COMM_WORLD, "** read time series %3.2e time %3.2e\n",
87  tit->first, tit->second);
88  }
89  }
90  if (r != 0) {
91  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
92  "*** ERROR file close unsuccessful");
93  }
94  readFile = 1;
96  }
97 
98  MoFEMErrorCode getForceScale(const double ts_t, double &scale) {
100  if (!fLg) {
101  scale = 1; // not scale at all, no history file
103  }
104  if (readFile == 0) {
105  SETERRQ(PETSC_COMM_SELF, 1, "data file not read");
106  }
107  scale = 0;
108  double t0 = 0, t1, s0 = tSeries[0], s1, dt;
109  std::map<double, double>::iterator tit = tSeries.begin();
110  for (; tit != tSeries.end(); tit++) {
111  if (tit->first > ts_t) {
112  t1 = tit->first;
113  s1 = tit->second;
114  dt = ts_t - t0;
115  scale = s0 + ((s1 - s0) / (t1 - t0)) * dt;
116  break;
117  }
118  t0 = tit->first;
119  s0 = tit->second;
120  scale = s0;
121  }
123  }
124 
125  /**
126  * @brief Scale force the right hand vector
127  *
128  * @param fe
129  * @param Nf
130  * @return MoFEMErrorCode
131  */
134  double scale;
135  const double ts_t = fe->ts_t;
136  CHKERR getForceScale(ts_t, scale);
137  Nf *= scale;
139  }
140 };
141 
143 
144  std::map<double, VectorDouble> tSeries;
146  string nAme;
147 
148  TimeAccelerogram(string name = "-my_accelerogram")
149  : readFile(0), debug(1), nAme(name) {
150 
151  ierr = timeData();
152  CHKERRABORT(PETSC_COMM_WORLD, ierr);
153  }
154 
157  char time_file_name[255];
158  PetscBool flg = PETSC_TRUE;
159  ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, nAme.c_str(),
160  time_file_name, 255, &flg);
161  CHKERRG(ierr);
162  if (flg != PETSC_TRUE) {
163  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
164  "*** ERROR %s (time_data FILE NEEDED)", nAme.c_str());
165  }
166  FILE *time_data = fopen(time_file_name, "r");
167  if (time_data == NULL) {
168  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
169  "*** ERROR data file < %s > open unsuccessful", time_file_name);
170  }
171  double no1 = 0.0;
172  VectorDouble no2(3);
173  tSeries[no1] = no2;
174  while (!feof(time_data)) {
175  int n =
176  fscanf(time_data, "%lf %lf %lf %lf", &no1, &no2[0], &no2[1], &no2[2]);
177  if (n < 0) {
178  fgetc(time_data);
179  continue;
180  }
181  if (n != 4) {
182  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
183  "*** ERROR read data file error (check input time data file) "
184  "{ n = %d }",
185  n);
186  }
187  tSeries[no1] = no2;
188  }
189  int r = fclose(time_data);
190  if (debug) {
191  std::map<double, VectorDouble>::iterator tit = tSeries.begin();
192  for (; tit != tSeries.end(); tit++) {
193  PetscPrintf(PETSC_COMM_WORLD,
194  "** read accelerogram %3.2e time %3.2e %3.2e %3.2e\n",
195  tit->first, tit->second[0], tit->second[1], tit->second[2]);
196  }
197  }
198  if (r != 0) {
199  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR file close unsuccessful");
200  }
201  readFile = 1;
203  }
204 
207  if (readFile == 0) {
208  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data file not read");
209  }
210  double ts_t = fe->ts_t;
211  VectorDouble acc(3);
212  VectorDouble acc0 = tSeries[0], acc1(3);
213  double t0 = 0, t1, dt;
214  std::map<double, VectorDouble>::iterator tit = tSeries.begin();
215  for (; tit != tSeries.end(); tit++) {
216  if (tit->first > ts_t) {
217  t1 = tit->first;
218  acc1 = tit->second;
219  dt = ts_t - t0;
220  acc = acc0 + ((acc1 - acc0) / (t1 - t0)) * dt;
221  break;
222  }
223  t0 = tit->first;
224  acc0 = tit->second;
225  acc = acc0;
226  }
227  Nf += acc;
229  }
230 };
231 
232 #endif // __TIMEFORCESCALE_HPP__
MoFEMErrorCode getForceScale(const double ts_t, double &scale)
TimeAccelerogram(string name="-my_accelerogram")
Force scale operator for reading two columns.
std::map< double, VectorDouble > tSeries
Class used to scale loads, f.e. in arc-length control.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:549
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &Nf)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
MoFEMErrorCode timeData()
#define CHKERR
Inline error check.
Definition: definitions.h:601
MoFEMErrorCode scaleNf(const FEMethod *fe, VectorDouble &Nf)
Scale force the right hand vector.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
std::map< double, double > tSeries
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
TimeForceScale(string name="-my_time_data_file", bool error_if_file_not_given=true)
MoFEMErrorCode timeData()