Skip to content

Commit 6d23972

Browse files
jfriedleinjfriedlein
committed
major update, change to standard S22
1 parent c43a034 commit 6d23972

File tree

3 files changed

+544
-931
lines changed

3 files changed

+544
-931
lines changed

‎DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh.h

+198-175
Original file line numberDiff line numberDiff line change
@@ -15,198 +15,221 @@
1515
// also contains enumerators as part of "enums::"
1616
#include "../numEx-helper_fnc.h"
1717

18-
using namespace dealii;
18+
#include "../numEx-baseClass.h"
1919

20-
namespace DENP_Laura_Abaqus_mesh
21-
/*
22-
* xxx
23-
*
24-
* CERTIFIED TO STANDARD xxx
20+
/**
21+
* @brief
22+
* CERTIFIED TO STANDARD S22
23+
*
24+
* @tparam dim
2525
*/
26+
template<int dim>
27+
class numEx_DENP_Laura_Abaqus_mesh : public numEx_class<dim>
2628
{
27-
// Name of the numerical example
28-
std::string numEx_name = "DENP_Laura_Abaqus_mesh";
29-
30-
// The loading direction: \n
31-
// In which coordinate direction the load shall be applied, so x/y/z.
32-
const unsigned int loading_direction = enums::y;
33-
34-
// The loaded faces:
35-
const enums::enum_boundary_ids id_boundary_load = enums::id_boundary_yPlus;
36-
const enums::enum_boundary_ids id_boundary_secondaryLoad = enums::id_boundary_xPlus;
37-
38-
// Characteristic body dimensions
39-
std::vector<double> body_dimensions (5);
40-
41-
// Some internal parameters
42-
struct parameterCollection
43-
{
44-
const double search_tolerance = 1e-12;
29+
public:
30+
std::string numEx_name() {
31+
return "DENP_Laura_Abaqus_mesh";
32+
};
33+
34+
unsigned int loading_direction() {
35+
return enums::y;
36+
};
37+
38+
std::vector< enums::enum_boundary_ids > id_boundary_loads()
39+
{
40+
std::vector< enums::enum_boundary_ids > id_boundary_loads_list (2);
41+
// The loaded faces:
42+
id_boundary_loads_list[enums::id_primary_load] = enums::id_boundary_yPlus;
43+
id_boundary_loads_list[enums::id_secondary_load] = enums::id_boundary_xPlus;
4544

45+
return id_boundary_loads_list;
46+
};
47+
48+
std::vector< types::manifold_id > manifold_ids()
49+
{
4650
const types::manifold_id manifold_id_right_radius = 11;
4751
const types::manifold_id manifold_id_left_radius = 10;
48-
};
49-
50-
const enums::enum_BC BC_xMinus = enums::BC_none;
51-
const enums::enum_BC BC_yPlus = enums::BC_x0_z0;
52-
const bool constrain_sideways_sliding_of_loaded_face = false;
53-
const enums::enum_BC BC_yMinus = enums::BC_fix;
54-
const enums::enum_BC BC_zMinus = enums::BC_sym;
55-
const enums::enum_notch_type notch_type = enums::notch_round;
56-
const bool notch_twice = true;
57-
const bool DENP_Laura = true;
58-
const bool DENP_Hagen = false;
59-
const bool SheStrip = false;
60-
61-
// Evaluation points: \n
62-
// We want points, one for the contraction of the center
63-
// and one for the contraction of the top face.
64-
// We don't know the coordinates yet, because the mesh has not yet been created.
65-
// So we fill the data in make_grid.
66-
// @todo We need \a dim here instead of 3, but dim is unkown at this place -> redesign
67-
std::vector< numEx::EvalPointClass<3> > eval_points_list (2, numEx::EvalPointClass<3>() );
68-
69-
/**
70-
* Apply the boundary conditions (support and load) on the given AffineConstraints \a constraints. \n
71-
* For the HyperRectangle that are three symmetry constraints on each plane (x=0, y=0, z=0) and the load on the \a id_boundary_load (for Dirichlet).
72-
*/
73-
template<int dim>
74-
void make_constraints ( AffineConstraints<double> &constraints, const FESystem<dim> &fe, DoFHandler<dim> &dof_handler_ref,
75-
const bool &apply_dirichlet_bc, double &current_load_increment,
76-
const Parameter::GeneralParameters &parameter )
52+
53+
return { manifold_id_left_radius, manifold_id_right_radius };
54+
}
55+
56+
void make_grid ( /*input-> */ const Parameter::GeneralParameters &parameter,
57+
/*output->*/ Triangulation<dim> &triangulation,
58+
std::vector<double> &body_dimensions,
59+
std::vector< numEx::EvalPointClass<3> > &eval_points_list,
60+
const std::string relativePath
61+
);
62+
63+
void make_constraints ( AffineConstraints<double> &constraints, const FESystem<dim> &fe, DoFHandler<dim> &dof_handler_ref,
64+
const bool &apply_dirichlet_bc, double &current_load_increment, const Parameter::GeneralParameters &parameter );
65+
66+
};
67+
68+
/**
69+
* Apply the boundary conditions (support and load) on the given AffineConstraints \a constraints. \n
70+
* For the HyperRectangle that are three symmetry constraints on each plane (x=0, y=0, z=0) and the load on the \a id_boundary_load (for Dirichlet).
71+
*/
72+
template<int dim>
73+
void numEx_DENP_Laura_Abaqus_mesh<dim>::make_constraints
74+
(
75+
AffineConstraints<double> &constraints, const FESystem<dim> &fe, DoFHandler<dim> &dof_handler_ref,
76+
const bool &apply_dirichlet_bc, double &current_load_increment,
77+
const Parameter::GeneralParameters &parameter
78+
)
79+
{
80+
const enums::enum_BC BC_xMinus = enums::BC_none;
81+
const enums::enum_BC BC_yPlus = enums::BC_x0_z0;
82+
const bool constrain_sideways_sliding_of_loaded_face = false;
83+
const enums::enum_BC BC_yMinus = enums::BC_fix;
84+
const enums::enum_BC BC_zMinus = enums::BC_sym;
85+
86+
// BC on x0 plane
87+
if ( BC_xMinus==enums::BC_sym )
88+
numEx::BC_apply( enums::id_boundary_xMinus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
89+
else if ( BC_xMinus==enums::BC_fix )
90+
numEx::BC_apply_fix( enums::id_boundary_xMinus, dof_handler_ref, fe, constraints );
91+
92+
// BC on y0 plane
93+
if ( BC_yMinus==enums::BC_fix )
7794
{
78-
parameterCollection parameters_internal;
79-
80-
// BC on x0 plane
81-
if ( BC_xMinus==enums::BC_sym )
82-
numEx::BC_apply( enums::id_boundary_xMinus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
83-
else if ( BC_xMinus==enums::BC_fix )
84-
numEx::BC_apply_fix( enums::id_boundary_xMinus, dof_handler_ref, fe, constraints );
85-
86-
// BC on y0 plane
87-
if ( BC_yMinus==enums::BC_fix )
88-
{
89-
// For compression we fix/clamp the Y0 plane, so it does not run away
90-
numEx::BC_apply_fix( enums::id_boundary_yMinus, dof_handler_ref, fe, constraints );
91-
}
92-
else if ( BC_yMinus==enums::BC_sym)
93-
numEx::BC_apply( enums::id_boundary_yMinus, enums::y, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
94-
95-
if ( BC_yPlus==enums::BC_x0_z0 )
96-
{
97-
numEx::BC_apply( enums::id_boundary_yPlus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
98-
numEx::BC_apply( enums::id_boundary_yPlus, enums::z, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
99-
}
100-
else if ( BC_yPlus==enums::BC_fix)
101-
numEx::BC_apply_fix( enums::id_boundary_yPlus, dof_handler_ref, fe, constraints );
102-
103-
// BC on z0 plane ...
104-
if ( dim==3 ) // ... only for 3D
105-
{
106-
// For compression we don't fix anything in the third direction, because y0 was already clamped.
107-
// @todo However, what about the upper part?
108-
if ( BC_zMinus == enums::BC_sym )
109-
numEx::BC_apply( enums::id_boundary_zMinus, enums::z, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
110-
}
111-
112-
// BC for the yPlus
113-
if ( constrain_sideways_sliding_of_loaded_face && BC_yPlus==enums::BC_x0 )
114-
numEx::BC_apply( enums::id_boundary_yPlus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
115-
116-
// BC for the load ...
117-
if ( parameter.driver == enums::Dirichlet ) // ... as Dirichlet only for Dirichlet as driver
118-
numEx::BC_apply( id_boundary_load, loading_direction, current_load_increment, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
95+
// For compression we fix/clamp the Y0 plane, so it does not run away
96+
numEx::BC_apply_fix( enums::id_boundary_yMinus, dof_handler_ref, fe, constraints );
11997
}
98+
else if ( BC_yMinus==enums::BC_sym)
99+
numEx::BC_apply( enums::id_boundary_yMinus, enums::y, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
100+
101+
if ( BC_yPlus==enums::BC_x0_z0 )
102+
{
103+
numEx::BC_apply( enums::id_boundary_yPlus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
104+
numEx::BC_apply( enums::id_boundary_yPlus, enums::z, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
105+
}
106+
else if ( BC_yPlus==enums::BC_fix)
107+
numEx::BC_apply_fix( enums::id_boundary_yPlus, dof_handler_ref, fe, constraints );
120108

121-
// 2D grid
122-
template <int dim>
123-
void make_grid( Triangulation<2> &triangulation, const Parameter::GeneralParameters &parameter )
109+
// BC on z0 plane ...
110+
if ( dim==3 ) // ... only for 3D
124111
{
112+
// For compression we don't fix anything in the third direction, because y0 was already clamped.
113+
// @todo However, what about the upper part?
114+
if ( BC_zMinus == enums::BC_sym )
115+
numEx::BC_apply( enums::id_boundary_zMinus, enums::z, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
116+
}
117+
118+
// BC for the yPlus
119+
if ( constrain_sideways_sliding_of_loaded_face && BC_yPlus==enums::BC_x0 )
120+
numEx::BC_apply( enums::id_boundary_yPlus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
125121

122+
// BC for the load ...
123+
if ( parameter.driver == enums::Dirichlet ) // ... as Dirichlet only for Dirichlet as driver
124+
{
125+
const std::vector< enums::enum_boundary_ids > id_boundary_loads_list = id_boundary_loads();
126+
numEx::BC_apply( id_boundary_loads_list[enums::id_primary_load], loading_direction(), current_load_increment, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
126127
}
128+
}
129+
130+
131+
// 2D grid
132+
// template <int dim>
133+
// void numEx_DENP_Laura_Abaqus_mesh<dim>::make_grid
134+
// (
135+
// /*input-> */const Parameter::GeneralParameters &parameter,
136+
// /*output->*/ Triangulation<dim> &triangulation,
137+
// std::vector<double> &body_dimensions,
138+
// std::vector< numEx::EvalPointClass<3> > &eval_points_list
139+
// )
140+
// {
141+
// AssertThrow(false, ExcMessage( numEx_name(), "make_grid<< Not yet available for 2D."));
142+
// }
127143

128-
// 3D grid
129-
template <int dim>
130-
void make_grid( Triangulation<3> &triangulation, const Parameter::GeneralParameters &parameter )
144+
145+
// 3D grid
146+
template <int dim>
147+
void numEx_DENP_Laura_Abaqus_mesh<dim>::make_grid
148+
(
149+
/*input-> */ const Parameter::GeneralParameters &parameter,
150+
/*output->*/ Triangulation<dim> &triangulation,
151+
std::vector<double> &body_dimensions,
152+
std::vector< numEx::EvalPointClass<3> > &eval_points_list,
153+
const std::string relativePath
154+
)
155+
{
156+
AssertThrow(dim==3, ExcMessage(numEx_name()+" << not yet available for 2D"));
157+
158+
const double search_tolerance = numEx_class<dim>::search_tolerance;
159+
160+
// Assign the dimensions of mesh and store them as characteristic lengths
161+
// @todo Currently hardcoded
162+
const double width = 18;
163+
body_dimensions[enums::x] = width;
164+
const double length = 60;
165+
body_dimensions[enums::y] = length;
166+
const double thickness = 1.5;
167+
body_dimensions[enums::z] = thickness;
168+
169+
// Import the Abaqus meshes
170+
GridIn<dim> grid_in;
171+
grid_in.attach_triangulation(triangulation);
172+
std::string path_to_inp;
173+
switch( parameter.refine_special )
174+
{
175+
case 1:
176+
path_to_inp = relativePath+"../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m1_792el.inp";
177+
break;
178+
case 2:
179+
path_to_inp = relativePath+"../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m2_2316el.inp";
180+
break;
181+
case 3:
182+
path_to_inp = relativePath+ "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m3_3584el.inp";
183+
break;
184+
case 4:
185+
path_to_inp = relativePath+"../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m4_3808el.inp";
186+
break;
187+
case 5:
188+
path_to_inp = relativePath+"../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m5_6040el.inp";
189+
break;
190+
case 6:
191+
path_to_inp = relativePath+"../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m6_10972el.inp";
192+
break;
193+
default:
194+
AssertThrow(false,ExcMessage("DENP_Laura_Abaqus_mesh::make_grid<< Not available mesh chosen via refine_special or default value used."));
195+
break;
196+
}
197+
std::ifstream input_file(path_to_inp);
198+
grid_in.read_abaqus(input_file);
199+
200+
// Clear all existing boundary ID's
201+
numEx::clear_boundary_IDs( triangulation );
202+
203+
// Set boundary IDs
204+
for (typename Triangulation<dim>::active_cell_iterator
205+
cell = triangulation.begin_active();
206+
cell != triangulation.end(); ++cell)
131207
{
132-
parameterCollection parameters_internal;
133-
const double search_tolerance = parameters_internal.search_tolerance;
134-
135-
// Assign the dimensions of the hyper rectangle and store them as characteristic lengths
136-
const double width = 18;
137-
body_dimensions[enums::x] = width;
138-
const double length = 60;
139-
body_dimensions[enums::y] = length;
140-
const double thickness = 1.5;
141-
body_dimensions[enums::z] = thickness;
142-
143-
GridIn<dim> grid_in;
144-
grid_in.attach_triangulation(triangulation);
145-
std::string path_to_inp;
146-
switch( parameter.refine_special )
147-
{
148-
case 1:
149-
path_to_inp = "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m1_792el.inp";
150-
break;
151-
case 2:
152-
path_to_inp = "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m2_2316el.inp";
153-
break;
154-
case 3:
155-
path_to_inp = "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m3_3584el.inp";
156-
break;
157-
case 4:
158-
path_to_inp = "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m4_3808el.inp";
159-
break;
160-
case 5:
161-
path_to_inp = "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m5_6040el.inp";
162-
break;
163-
case 6:
164-
path_to_inp = "../Numerical_examples_in_dealii/DENP_Laura_Abaqus_mesh/DENP_Laura_Abaqus_mesh_m6_10972el.inp";
165-
break;
166-
default:
167-
AssertThrow(false,ExcMessage("DENP_Laura_Abaqus_mesh::make_grid<< Not available mesh chosen via refine_special or default value used."));
168-
break;
169-
}
170-
std::ifstream input_file(path_to_inp);
171-
172-
grid_in.read_abaqus(input_file);
173-
174-
// Clear all existing boundary ID's
175-
numEx::clear_boundary_IDs( triangulation );
176-
177-
// Set boundary IDs
178-
for (typename Triangulation<dim>::active_cell_iterator
179-
cell = triangulation.begin_active();
180-
cell != triangulation.end(); ++cell)
208+
for (unsigned int face=0; face < GeometryInfo<dim>::faces_per_cell; ++face)
181209
{
182-
for (unsigned int face=0; face < GeometryInfo<dim>::faces_per_cell; ++face)
210+
if (cell->face(face)->at_boundary())
183211
{
184-
if (cell->face(face)->at_boundary())
185-
{
186-
// Set boundary IDs
187-
if (std::abs(cell->face(face)->center()[enums::x] - 0.0) < search_tolerance)
188-
cell->face(face)->set_boundary_id(enums::id_boundary_xMinus);
189-
else if ( std::abs(cell->face(face)->center()[enums::x] - body_dimensions[enums::x] ) < search_tolerance)
190-
cell->face(face)->set_boundary_id(enums::id_boundary_xPlus);
191-
else if (std::abs(cell->face(face)->center()[enums::y] - 0.0) < search_tolerance)
192-
cell->face(face)->set_boundary_id(enums::id_boundary_yMinus);
193-
else if (std::abs(cell->face(face)->center()[enums::y] - body_dimensions[enums::y]) < search_tolerance)
194-
cell->face(face)->set_boundary_id(enums::id_boundary_yPlus);
195-
else if ( std::abs(cell->face(face)->center()[enums::z] - 0.0) < search_tolerance)
196-
cell->face(face)->set_boundary_id(enums::id_boundary_zMinus);
197-
else if ( std::abs(cell->face(face)->center()[enums::z] - body_dimensions[enums::z]) < search_tolerance)
198-
cell->face(face)->set_boundary_id(enums::id_boundary_zPlus);
199-
}
212+
// Set boundary IDs
213+
if (std::abs(cell->face(face)->center()[enums::x] - 0.0) < search_tolerance)
214+
cell->face(face)->set_boundary_id(enums::id_boundary_xMinus);
215+
else if ( std::abs(cell->face(face)->center()[enums::x] - body_dimensions[enums::x] ) < search_tolerance)
216+
cell->face(face)->set_boundary_id(enums::id_boundary_xPlus);
217+
else if (std::abs(cell->face(face)->center()[enums::y] - 0.0) < search_tolerance)
218+
cell->face(face)->set_boundary_id(enums::id_boundary_yMinus);
219+
else if (std::abs(cell->face(face)->center()[enums::y] - body_dimensions[enums::y]) < search_tolerance)
220+
cell->face(face)->set_boundary_id(enums::id_boundary_yPlus);
221+
else if ( std::abs(cell->face(face)->center()[enums::z] - 0.0) < search_tolerance)
222+
cell->face(face)->set_boundary_id(enums::id_boundary_zMinus);
223+
else if ( std::abs(cell->face(face)->center()[enums::z] - body_dimensions[enums::z]) < search_tolerance)
224+
cell->face(face)->set_boundary_id(enums::id_boundary_zPlus);
200225
}
201-
} // end for(cell)
226+
}
227+
} // end for(cell)
202228

203-
// Evaluation points and the related list of them
204-
numEx::EvalPointClass<3> eval_center ( Point<3>(width-2.5,0,0), enums::x );
205-
numEx::EvalPointClass<3> eval_top ( Point<3>(body_dimensions[enums::x],body_dimensions[enums::y],0), enums::x );
229+
// Evaluation points
230+
const numEx::EvalPointClass<3> eval_top ( Point<3>(body_dimensions[enums::x],body_dimensions[enums::y],0), enums::y );
231+
eval_points_list[enums::eval_point_0] = eval_top;
206232

207-
eval_points_list = {eval_center,eval_top};
208-
209-
// Output the triangulation as eps or inp
210-
// numEx::output_triangulation( triangulation, enums::output_eps_as_2D, numEx_name, true );
211-
}
233+
// Output the triangulation as eps or inp
234+
// numEx::output_triangulation( triangulation, enums::output_eps_as_2D, numEx_name, true );
212235
}

0 commit comments

Comments
 (0)