-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathHyperCube.h
174 lines (148 loc) · 6.87 KB
/
HyperCube.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
// Numerical example helper function (required, can be downloaded from https://github.com/jfriedlein/Numerical_examples_in_dealii)
// also contains enumerators as part of "enums::"
#include "numEx-helper_fnc.h"
#include "numEx-baseClass.h"
/**
* @brief
* CERTIFIED TO STANDARD S22
*
* @tparam dim
*/
template<int dim>
class numEx_HyperCube : public numEx_class<dim>
{
public:
std::string numEx_name() {
return "HyperCube";
};
unsigned int loading_direction() {
return enums::y;
};
std::vector< enums::enum_boundary_ids > id_boundary_loads()
{
// The loaded faces:
const enums::enum_boundary_ids id_boundary_load = enums::id_boundary_yPlus;
const enums::enum_boundary_ids id_boundary_secondaryLoad = enums::id_boundary_xPlus;
std::vector< enums::enum_boundary_ids > id_boundary_loads_list (2);
id_boundary_loads_list[enums::id_primary_load] = id_boundary_load;
id_boundary_loads_list[enums::id_secondary_load] = id_boundary_secondaryLoad;
return id_boundary_loads_list;
};
void make_grid ( /*input-> */ const Parameter::GeneralParameters ¶meter,
/*output->*/ Triangulation<dim> &triangulation,
std::vector<double> &body_dimensions,
std::vector< numEx::EvalPointClass<3> > &eval_points_list,
const std::string relativePath
);
void make_constraints ( AffineConstraints<double> &constraints, const FESystem<dim> &fe, DoFHandler<dim> &dof_handler_ref,
const bool &apply_dirichlet_bc, double ¤t_load_increment, const Parameter::GeneralParameters ¶meter );
};
// HyperCube grid: 2D and 3D
template<int dim>
void numEx_HyperCube<dim>::make_grid
(
/*input-> */ const Parameter::GeneralParameters ¶meter,
/*output->*/ Triangulation<dim> &triangulation,
std::vector<double> &body_dimensions,
std::vector< numEx::EvalPointClass<3> > &eval_points_list,
const std::string relativePath
)
{
const double search_tolerance = numEx_class<dim>::search_tolerance;
const double width = parameter.width;
// Assign the characteristic dimensions of the cube
body_dimensions[enums::x] = width;
body_dimensions[enums::y] = width;
body_dimensions[enums::z] = width;
// Evaluation points and the related list of them
const numEx::EvalPointClass<3> eval_top ( Point<3>(0,body_dimensions[enums::y],0), enums::y );
eval_points_list[enums::eval_point_0] = eval_top;
const numEx::EvalPointClass<3> eval_path_start ( Point<3>(0,0,0), enums::coord_none );
const numEx::EvalPointClass<3> eval_path_end = eval_top;
eval_points_list[enums::eval_path0_start] = eval_path_start;
eval_points_list[enums::eval_path0_end ] = eval_path_end;
// Create the triangulation
GridGenerator::hyper_cube(triangulation,0,width);
// Clear all existing boundary ID's
numEx::clear_boundary_IDs( triangulation );
// Set the boundary IDs
for ( typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell )
{
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary())
{
if (std::abs(cell->face(face)->center()[enums::x] - 0.0) < search_tolerance)
{
cell->face(face)->set_boundary_id(enums::id_boundary_xMinus);
}
else if (std::abs(cell->face(face)->center()[enums::x] - body_dimensions[enums::x]) < search_tolerance)
{
cell->face(face)->set_boundary_id(enums::id_boundary_xPlus);
}
else if (std::abs(cell->face(face)->center()[enums::y] - 0.0) < search_tolerance)
{
cell->face(face)->set_boundary_id(enums::id_boundary_yMinus);
}
else if (std::abs(cell->face(face)->center()[enums::y] - body_dimensions[enums::y]) < search_tolerance)
{
cell->face(face)->set_boundary_id(enums::id_boundary_yPlus);
}
else if (dim==3 && std::abs(cell->face(face)->center()[enums::z] - 0.0) < search_tolerance)
{
cell->face(face)->set_boundary_id(enums::id_boundary_zMinus);
}
else if (dim==3 && std::abs(cell->face(face)->center()[enums::z] - body_dimensions[enums::z]) < search_tolerance)
{
cell->face(face)->set_boundary_id(enums::id_boundary_zPlus);
}
else
{
// There are only 6 faces for a cube in 3D, so if we missed one, something went terribly wrong
AssertThrow(false, ExcMessage( numEx_name()+" - make_grid<< Found an unidentified face at the boundary. "
"Maybe it slipt through the assignment or that face is simply not needed. "
"So either check the implementation or comment this line in the code") );
// Just to get rid of the unused variable warning
std::cout << "Relative path " << relativePath << std::endl;
}
}
}
// Refinement
triangulation.refine_global( parameter.nbr_global_refinements );
// Output the triangulation as eps or inp
//numEx::output_triangulation( triangulation, enums::output_eps, numEx_name );
}
/**
* Apply the boundary conditions (support and load) on the given AffineConstraints \a constraints. \n
* For the HyperCube 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).
* Alternatively, we apply the rigid body motion for contact.
* @todo Change setup (see HyperR) where we define boundary for BC_xMinus and then use if ( BC_xMinus==... ), use the value of BC_xMinus directly
*/
template<int dim>
void numEx_HyperCube<dim>::make_constraints ( AffineConstraints<double> &constraints, const FESystem<dim> &fe, DoFHandler<dim> &dof_handler_ref,
const bool &apply_dirichlet_bc, double ¤t_load_increment, const Parameter::GeneralParameters ¶meter )
{
// BC on x0 plane
numEx::BC_apply( enums::id_boundary_xMinus, enums::x, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
// BC on y0 plane
numEx::BC_apply( enums::id_boundary_yMinus, enums::y, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
// Fix the bottom face
//numEx::BC_apply_fix( enums::id_boundary_yMinus, dof_handler_ref, fe, constraints );
// BC on z0 plane ...
if ( dim==3 ) // ... only for 3D
numEx::BC_apply( enums::id_boundary_zMinus, enums::z, 0, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
// BC for the load ...
if ( parameter.driver == enums::Dirichlet ) // ... as Dirichlet only for Dirichlet as driver, alternatively ...
{
const std::vector< enums::enum_boundary_ids > id_boundary_loads_list = id_boundary_loads();
const unsigned int loading_dir = loading_direction();
numEx::BC_apply( id_boundary_loads_list[enums::id_primary_load], loading_dir, current_load_increment, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
//numEx::BC_apply( id_boundary_loads_list[1], enums::x, current_load_increment, apply_dirichlet_bc, dof_handler_ref, fe, constraints );
}
else if ( parameter.driver == enums::Contact ) // ... as contact
{
//if (apply_dirichlet_bc == true )
// rigid_wall->move( current_load_increment );
}
}