13 #if !defined(KRATOS_INITIAL_STRESS_2D_UTILITIES )
14 #define KRATOS_INITIAL_STRESS_2D_UTILITIES
69 this->ComputeCellMatrixDimensions(AuxVariables,rCurrentModelPart);
78 std::fstream initial_stresses_mdpa;
79 std::string initial_stress_mdpa_path = rParameters[
"initial_input_filename"].
GetString() +
".mdpa";
81 initial_stresses_mdpa.open(initial_stress_mdpa_path.c_str(),
std::fstream::out | std::fstream::trunc);
82 initial_stresses_mdpa.precision(15);
84 initial_stresses_mdpa <<
"Begin Properties 0" << std::endl;
85 initial_stresses_mdpa <<
"End Properties" << std::endl << std::endl;
87 const int NNodes =
static_cast<int>(rCurrentModelPart.
Nodes().size());
88 ModelPart::NodesContainerType::iterator node_begin = rCurrentModelPart.
NodesBegin();
89 initial_stresses_mdpa <<
"Begin Nodes" << std::endl;
91 for(
int i = 0;
i < NNodes;
i++) {
92 ModelPart::NodesContainerType::iterator it_node = node_begin +
i;
93 initial_stresses_mdpa << it_node->Id() <<
" " <<
94 it_node->X0() <<
" " <<
95 it_node->Y0() <<
" " <<
96 it_node->Z0() << std::endl;
98 initial_stresses_mdpa <<
"End Nodes" << std::endl << std::endl;
100 this->WriteLinearElements(initial_stresses_mdpa,rCurrentModelPart);
102 this->CalculateNodalStresses(rCurrentModelPart);
103 Matrix InitialStressTensor(3,3);
104 initial_stresses_mdpa <<
"Begin NodalData INITIAL_STRESS_TENSOR" << std::endl;
106 for(
int i = 0;
i < NNodes;
i++) {
107 ModelPart::NodesContainerType::iterator it_node = node_begin +
i;
108 noalias(InitialStressTensor) = it_node->FastGetSolutionStepValue(NODAL_EFFECTIVE_STRESS_TENSOR);
109 initial_stresses_mdpa << it_node->Id() <<
" 0 [2,2] ((" <<
110 InitialStressTensor(0,0) <<
"," << InitialStressTensor(0,1) <<
"),(" <<
111 InitialStressTensor(1,0) <<
"," << InitialStressTensor(1,1) <<
"))" << std::endl;
113 initial_stresses_mdpa <<
"End NodalData" << std::endl;
115 initial_stresses_mdpa.close();
130 std::vector< std::vector< std::vector<Element::Pointer> > > ElementOldCellMatrix;
131 ElementOldCellMatrix.resize(AuxVariables.
NRows);
132 for(
int i = 0;
i < AuxVariables.
NRows;
i++) ElementOldCellMatrix[
i].resize(AuxVariables.
NColumns);
139 int NElems =
static_cast<int>(rInitialModelPart.
Elements().size());
140 ModelPart::ElementsContainerType::iterator el_begin = rInitialModelPart.
ElementsBegin();
142 #pragma omp parallel for private(X_me,Y_me,PointsNumber)
143 for(
int i = 0;
i < NElems;
i++)
145 ModelPart::ElementsContainerType::iterator itElemOld = el_begin +
i;
147 double X_left = itElemOld->GetGeometry().GetPoint(0).X0();
148 double X_right = X_left;
149 double Y_top = itElemOld->GetGeometry().GetPoint(0).Y0();
150 double Y_bot = Y_top;
151 PointsNumber = itElemOld->GetGeometry().PointsNumber();
153 for(
int j = 1;
j < PointsNumber;
j++)
155 X_me = itElemOld->GetGeometry().GetPoint(
j).X0();
156 Y_me = itElemOld->GetGeometry().GetPoint(
j).Y0();
158 if(X_me > X_right) X_right = X_me;
159 else if(X_me < X_left) X_left = X_me;
160 if(Y_me > Y_top) Y_top = Y_me;
161 else if(Y_me < Y_bot) Y_bot = Y_me;
169 if(Column_left < 0) Column_left = 0;
170 else if(Column_left >= AuxVariables.
NColumns) Column_left = AuxVariables.
NColumns-1;
171 if(Column_right >= AuxVariables.
NColumns) Column_right = AuxVariables.
NColumns-1;
172 else if(Column_right < 0) Column_right = 0;
174 if(Row_top < 0) Row_top = 0;
175 else if(Row_top >= AuxVariables.
NRows) Row_top = AuxVariables.
NRows-1;
176 if(Row_bot >= AuxVariables.
NRows) Row_bot = AuxVariables.
NRows-1;
177 else if(Row_bot < 0) Row_bot = 0;
179 for(
int k = Row_top;
k <= Row_bot;
k++)
181 for(
int l = Column_left;
l <= Column_right;
l++)
185 ElementOldCellMatrix[
k][
l].push_back((*(itElemOld.base())));
192 const int NNodes =
static_cast<int>(rCurrentModelPart.
Nodes().size());
193 ModelPart::NodesContainerType::iterator node_begin = rCurrentModelPart.
NodesBegin();
198 #pragma omp parallel for private(X_me,Y_me,PointsNumber,GlobalCoordinates,LocalCoordinates)
199 for(
int i = 0;
i < NNodes;
i++)
201 ModelPart::NodesContainerType::iterator itNodeNew = node_begin +
i;
203 X_me = itNodeNew->X0();
204 Y_me = itNodeNew->Y0();
210 else if(Column < 0) Column = 0;
211 if(Row >= AuxVariables.
NRows) Row = AuxVariables.
NRows-1;
212 else if(Row < 0) Row = 0;
214 noalias(GlobalCoordinates) = itNodeNew->Coordinates();
215 bool IsInside =
false;
216 Element::Pointer pElementOld;
218 for(
unsigned int m = 0;
m < (ElementOldCellMatrix[Row][Column]).size();
m++)
220 pElementOld = ElementOldCellMatrix[Row][Column][
m];
221 IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates);
226 for(
unsigned int m = 0;
m < (ElementOldCellMatrix[Row][Column]).size();
m++)
228 pElementOld = ElementOldCellMatrix[Row][Column][
m];
229 IsInside = pElementOld->GetGeometry().IsInside(GlobalCoordinates,LocalCoordinates,1.0e-5);
233 if(IsInside ==
false)
234 std::cout <<
"ERROR!!, NONE OF THE OLD ELEMENTS CONTAINS NODE: " << itNodeNew->Id() << std::endl;
236 PointsNumber = pElementOld->GetGeometry().PointsNumber();
238 Vector ShapeFunctionsValuesVector(PointsNumber);
239 pElementOld->GetGeometry().ShapeFunctionsValues(ShapeFunctionsValuesVector,LocalCoordinates);
242 std::vector<Vector> ComponentsNodalVariableVector(3);
243 for(
int j = 0;
j < 3;
j++) {
244 ComponentsNodalVariableVector[
j].resize(PointsNumber);
246 Vector nodal_initial_stress_vector(3);
247 Matrix nodal_initial_stress_tensor(2,2);
248 for(
int j = 0;
j < PointsNumber;
j++) {
249 noalias(nodal_initial_stress_tensor) = pElementOld->GetGeometry().GetPoint(
j).FastGetSolutionStepValue(INITIAL_STRESS_TENSOR);
251 for (
int k= 0;
k < 3;
k++) {
252 ComponentsNodalVariableVector[
k][
j] = nodal_initial_stress_vector[
k];
255 for (
int k= 0;
k < 3;
k++) {
256 nodal_initial_stress_vector[
k] =
inner_prod(ShapeFunctionsValuesVector,ComponentsNodalVariableVector[
k]);
259 Matrix& rNodalStress = itNodeNew->FastGetSolutionStepValue(INITIAL_STRESS_TENSOR);
260 if(rNodalStress.size1() != 2)
261 rNodalStress.
resize(2,2,
false);
262 noalias(rNodalStress) = nodal_initial_stress_tensor;
272 void ComputeCellMatrixDimensions(
273 UtilityVariables& rAuxVariables,
278 std::vector<double> X_max_partition(NumThreads);
279 std::vector<double> X_min_partition(NumThreads);
280 std::vector<double> Y_max_partition(NumThreads);
281 std::vector<double> Y_min_partition(NumThreads);
283 const int NNodes =
static_cast<int>(rModelPart.
Nodes().size());
284 ModelPart::NodesContainerType::iterator node_begin = rModelPart.
NodesBegin();
290 X_max_partition[
k] = node_begin->X();
291 X_min_partition[
k] = X_max_partition[
k];
292 Y_max_partition[
k] = node_begin->Y();
293 Y_min_partition[
k] = Y_max_partition[
k];
298 for(
int i = 0;
i < NNodes;
i++)
300 ModelPart::NodesContainerType::iterator itNode = node_begin +
i;
305 if( X_me > X_max_partition[
k] ) X_max_partition[
k] = X_me;
306 else if( X_me < X_min_partition[
k] ) X_min_partition[
k] = X_me;
308 if( Y_me > Y_max_partition[
k] ) Y_max_partition[
k] = Y_me;
309 else if( Y_me < Y_min_partition[
k] ) Y_min_partition[
k] = Y_me;
313 rAuxVariables.X_max = X_max_partition[0];
314 rAuxVariables.X_min = X_min_partition[0];
315 rAuxVariables.Y_max = Y_max_partition[0];
316 rAuxVariables.Y_min = Y_min_partition[0];
318 for(
unsigned int i=1;
i < NumThreads;
i++)
320 if(X_max_partition[
i] > rAuxVariables.X_max) rAuxVariables.X_max = X_max_partition[
i];
321 if(X_min_partition[
i] < rAuxVariables.X_min) rAuxVariables.X_min = X_min_partition[
i];
322 if(Y_max_partition[
i] > rAuxVariables.Y_max) rAuxVariables.Y_max = Y_max_partition[
i];
323 if(Y_min_partition[
i] < rAuxVariables.Y_min) rAuxVariables.Y_min = Y_min_partition[
i];
327 double AverageElementLength = 0.0;
329 int NElems =
static_cast<int>(rModelPart.
Elements().size());
330 ModelPart::ElementsContainerType::iterator el_begin = rModelPart.
ElementsBegin();
332 #pragma omp parallel for reduction(+:AverageElementLength)
333 for(
int i = 0;
i < NElems;
i++)
335 ModelPart::ElementsContainerType::iterator itElem = el_begin +
i;
337 AverageElementLength += itElem->GetGeometry().Length();
340 AverageElementLength = AverageElementLength/NElems;
343 rAuxVariables.NRows =
int((rAuxVariables.Y_max-rAuxVariables.Y_min)/(AverageElementLength));
344 rAuxVariables.NColumns =
int((rAuxVariables.X_max-rAuxVariables.X_min)/(AverageElementLength));
345 if(rAuxVariables.NRows <= 0) rAuxVariables.NRows = 1;
346 if(rAuxVariables.NColumns <= 0) rAuxVariables.NColumns = 1;
347 rAuxVariables.RowSize = (rAuxVariables.Y_max-rAuxVariables.Y_min)/rAuxVariables.NRows;
348 rAuxVariables.ColumnSize = (rAuxVariables.X_max-rAuxVariables.X_min)/rAuxVariables.NColumns;
353 void CalculateNodalStresses(ModelPart& rCurrentModelPart) {
356 const int NNodes =
static_cast<int>(rCurrentModelPart.Nodes().size());
357 ModelPart::NodesContainerType::iterator node_begin = rCurrentModelPart.NodesBegin();
359 #pragma omp parallel for
360 for(
int i = 0;
i < NNodes;
i++)
362 ModelPart::NodesContainerType::iterator it_node = node_begin +
i;
363 it_node->FastGetSolutionStepValue(NODAL_AREA) = 0.0;
364 Matrix& rNodalStress = it_node->FastGetSolutionStepValue(NODAL_EFFECTIVE_STRESS_TENSOR);
365 if(rNodalStress.size1() != 3)
366 rNodalStress.
resize(3,3,
false);
371 const ProcessInfo& r_current_process_info = rCurrentModelPart.GetProcessInfo();
372 const auto it_elem_begin = rCurrentModelPart.ElementsBegin();
373 #pragma omp parallel for
374 for(
int i=0; i<static_cast<int>(rCurrentModelPart.Elements().size()); ++
i) {
375 auto it_elem = it_elem_begin +
i;
376 it_elem->FinalizeSolutionStep(r_current_process_info);
380 #pragma omp parallel for
381 for(
int n = 0;
n < NNodes;
n++)
383 ModelPart::NodesContainerType::iterator it_node = node_begin +
n;
384 const double& NodalArea = it_node->FastGetSolutionStepValue(NODAL_AREA);
385 if (NodalArea>1.0e-20)
387 const double InvNodalArea = 1.0/NodalArea;
388 Matrix& rNodalStress = it_node->FastGetSolutionStepValue(NODAL_EFFECTIVE_STRESS_TENSOR);
389 for(
unsigned int i = 0;
i<3;
i++)
391 for(
unsigned int j = 0;
j<3;
j++)
393 rNodalStress(
i,
j) *= InvNodalArea;
402 void WriteLinearElements(std::fstream& rInitialStressesMdpa, ModelPart& rCurrentModelPart) {
404 int NElems =
static_cast<int>(rCurrentModelPart.Elements().size());
405 ModelPart::ElementsContainerType::iterator el_begin = rCurrentModelPart.ElementsBegin();
407 int PointsNumber = el_begin->GetGeometry().PointsNumber();
408 if(PointsNumber == 3){
409 rInitialStressesMdpa <<
"Begin Elements Element2D3N" << std::endl;
411 for(
int i = 0;
i < NElems;
i++)
413 ModelPart::ElementsContainerType::iterator it_elem = el_begin +
i;
414 rInitialStressesMdpa << it_elem->Id() <<
" 0 " <<
415 it_elem->GetGeometry().GetPoint(0).Id() <<
" " <<
416 it_elem->GetGeometry().GetPoint(1).Id() <<
" " <<
417 it_elem->GetGeometry().GetPoint(2).Id() << std::endl;
420 rInitialStressesMdpa <<
"Begin Elements Element2D4N" << std::endl;
422 for(
int i = 0;
i < NElems;
i++)
424 ModelPart::ElementsContainerType::iterator it_elem = el_begin +
i;
425 rInitialStressesMdpa << it_elem->Id() <<
" 0 " <<
426 it_elem->GetGeometry().GetPoint(0).Id() <<
" " <<
427 it_elem->GetGeometry().GetPoint(1).Id() <<
" " <<
428 it_elem->GetGeometry().GetPoint(2).Id() <<
" " <<
429 it_elem->GetGeometry().GetPoint(3).Id() << std::endl;
432 rInitialStressesMdpa <<
"End Elements" << std::endl << std::endl;
Definition: initial_stress_2D_utilities.hpp:37
virtual ~InitialStress2DUtilities()
Destructor.
Definition: initial_stress_2D_utilities.hpp:60
void NodalVariablesMapping(const UtilityVariables &AuxVariables, ModelPart &rInitialModelPart, ModelPart &rCurrentModelPart)
Member Variables.
Definition: initial_stress_2D_utilities.hpp:124
InitialStress2DUtilities()
Constructor.
Definition: initial_stress_2D_utilities.hpp:55
void TransferInitialStresses(ModelPart &rInitialModelPart, ModelPart &rCurrentModelPart)
Definition: initial_stress_2D_utilities.hpp:64
KRATOS_CLASS_POINTER_DEFINITION(InitialStress2DUtilities)
void SaveInitialStresses(Parameters &rParameters, ModelPart &rCurrentModelPart)
Definition: initial_stress_2D_utilities.hpp:76
void resize(std::size_t NewSize1, std::size_t NewSize2, bool preserve=0)
Definition: amatrix_interface.h:224
static TMatrixType StressVectorToTensor(const TVector &rStressVector)
Transforms a stess vector into a matrix. Stresses are assumed to be stored in the following way: for...
Definition: math_utils.h:1174
static TVector StressTensorToVector(const TMatrixType &rStressTensor, SizeType rSize=0)
Transforms a given symmetric Stress Tensor to Voigt Notation:
Definition: math_utils.h:1399
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
ElementIterator ElementsBegin(IndexType ThisIndex=0)
Definition: model_part.h:1169
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
ElementsContainerType & Elements(IndexType ThisIndex=0)
Definition: model_part.h:1189
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
This class provides to Kratos a data structure for I/O based on the standard of JSON.
Definition: kratos_parameters.h:59
std::string GetString() const
This method returns the string contained in the current Parameter.
Definition: kratos_parameters.cpp:684
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Internals::Matrix< double, AMatrix::dynamic, AMatrix::dynamic > Matrix
Definition: amatrix_interface.h:470
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
REACTION_CHECK_STIFFNESS_FACTOR int
Definition: contact_structural_mechanics_application_variables.h:75
ProcessInfo
Definition: edgebased_PureConvection.py:116
out
Definition: isotropic_damage_automatic_differentiation.py:200
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
int k
Definition: quadrature.py:595
int j
Definition: quadrature.py:648
int m
Definition: run_marine_rain_substepping.py:8
integer i
Definition: TensorModule.f:17
integer l
Definition: TensorModule.f:17
Basic Structs for the utility -----------------------------------------------------------------------...
Definition: initial_stress_2D_utilities.hpp:44
int NRows
Definition: initial_stress_2D_utilities.hpp:46
int NColumns
Definition: initial_stress_2D_utilities.hpp:46
double Y_max
Definition: initial_stress_2D_utilities.hpp:45
double ColumnSize
Definition: initial_stress_2D_utilities.hpp:47
double X_min
Definition: initial_stress_2D_utilities.hpp:45
double X_max
Definition: initial_stress_2D_utilities.hpp:45
double Y_min
Definition: initial_stress_2D_utilities.hpp:45
double RowSize
Definition: initial_stress_2D_utilities.hpp:47