CArl
Code Arlequin / C++ implementation
assemble_functions_elasticity_anisotropy_3D.cpp
Go to the documentation of this file.
1 /*
2  * assemble_functions_elasticity_anisotropy.cpp
3  *
4  * Created on: Sep 3, 2016
5  * Author: Thiago Milanetto Schlittler
6  */
7 
8 #include "assemble_functions_elasticity_anisotropy_3D.h"
9 
10 void Update_SubK_anisotropic( libMesh::DenseSubMatrix<libMesh::Number>& SubK,
11  unsigned int qp,
12  unsigned int C_i,
13  unsigned int C_k,
14  const std::vector<std::vector<libMesh::RealGradient> >& dphi,
15  const unsigned int n_components,
16  const unsigned int n_u_dofs,
17  const std::vector<libMesh::Real>& JxW,
18  carl::anisotropic_elasticity_tensor_cubic_sym& anisotropy_obj_input,
19  int grain_idx,
20  double cte
21  )
22 {
23  for (unsigned int iii=0; iii<n_u_dofs; iii++)
24  {
25  for (unsigned int jjj=0; jjj<n_u_dofs; jjj++)
26  {
27  for(unsigned int C_j=0; C_j<n_components; C_j++)
28  {
29  for(unsigned int C_l=0; C_l<n_components; C_l++)
30  {
31  SubK(iii,jjj) += cte * JxW[qp]*(anisotropy_obj_input.eval_elasticity_tensor(C_i,C_j,C_k,C_l,grain_idx) * dphi[iii][qp](C_j)*dphi[jjj][qp](C_l));
32  }
33  }
34  }
35  }
36 };
37 
38 void assemble_elasticity_anisotropic_with_weight( libMesh::EquationSystems& es,
39  const std::string& system_name,
40  carl::weight_parameter_function& weight_mask,
41  carl::anisotropic_elasticity_tensor_cubic_sym& anisotropy_obj_input)
42 {
43  libmesh_assert_equal_to (system_name, "Elasticity");
44 
45  libMesh::PerfLog perf_log ("Matrix Assembly (Heterogeneous elasticity)",MASTER_bPerfLog_assemble_fem);
46 
47  perf_log.push("Preamble");
48  // Set up mesh
49  const libMesh::MeshBase& mesh = es.get_mesh();
50 
51  const unsigned int dim = mesh.mesh_dimension();
52 
53  // - Set up physical properties system ------------------------------------
54  int localIdx;
55 
56  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
57  const unsigned int idx_var = physical_param_system.variable_number ("Index");
58 
59  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
60  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
61 
62  // - Set up elasticity system ---------------------------------------------
63  libMesh::LinearImplicitSystem& system = es.get_system<libMesh::LinearImplicitSystem>("Elasticity");
64 
65  const unsigned int n_components = 3;
66  const unsigned int u_var = system.variable_number ("u");
67  const unsigned int v_var = system.variable_number ("v");
68  const unsigned int w_var = system.variable_number ("w");
69 
70  const libMesh::DofMap& dof_map = system.get_dof_map();
71  libMesh::FEType fe_type = dof_map.variable_type(u_var);
72 
73  // Set up pointers to FEBase's of dimension dim and FE type fe_type
74  // -> 3D elements
75  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
76  libMesh::QGauss qrule (dim, fe_type.default_quadrature_order());
77  fe->attach_quadrature_rule (&qrule);
78 
79  // -> Faces
80  libMesh::UniquePtr<libMesh::FEBase> fe_face (libMesh::FEBase::build(dim, fe_type));
81  libMesh::QGauss qface(dim-1, fe_type.default_quadrature_order());
82  fe_face->attach_quadrature_rule (&qface);
83 
84  // Jacobian
85  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
86 
87  // Shape functions
88  const std::vector<std::vector<libMesh::Real> >& phi = fe->get_phi();
89 
90  // Shape functions derivatives
91  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
92 
93  // Quadrature points
94  const std::vector<libMesh::Point>& qp_points = fe->get_xyz();
95 
96  // Weights for the Arlequin method
97  double alpha_micro = 1;
98 
99  libMesh::DenseMatrix<libMesh::Number> Ke;
100  libMesh::DenseVector<libMesh::Number> Fe;
101 
102  libMesh::DenseSubMatrix<libMesh::Number>
103  Kuu(Ke), Kuv(Ke), Kuw(Ke),
104  Kvu(Ke), Kvv(Ke), Kvw(Ke),
105  Kwu(Ke), Kwv(Ke), Kww(Ke);
106 
107  libMesh::DenseSubVector<libMesh::Number>
108  Fu(Fe),
109  Fv(Fe),
110  Fw(Fe);
111 
112  std::vector<libMesh::dof_id_type> dof_indices;
113  std::vector<libMesh::dof_id_type> dof_indices_u;
114  std::vector<libMesh::dof_id_type> dof_indices_v;
115  std::vector<libMesh::dof_id_type> dof_indices_w;
116 
117  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
118  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
119 
120  perf_log.pop("Preamble");
121 
122  // For each element
123  for ( ; el != end_el; ++el)
124  {
125  // Get its pointer
126  const libMesh::Elem* elem = *el;
127 
128  perf_log.push("Define DoF");
129  // The total DoF indices, and those associated to each variable
130  dof_map.dof_indices (elem, dof_indices);
131  dof_map.dof_indices (elem, dof_indices_u, u_var);
132  dof_map.dof_indices (elem, dof_indices_v, v_var);
133  dof_map.dof_indices (elem, dof_indices_w, w_var);
134 
135  const unsigned int n_dofs = dof_indices.size();
136  const unsigned int n_u_dofs = dof_indices_u.size();
137  const unsigned int n_v_dofs = dof_indices_v.size();
138  const unsigned int n_w_dofs = dof_indices_w.size();
139 
140  perf_log.pop("Define DoF");
141 
142  perf_log.push("Define physical params");
143  // The DoF and values of the physical system
144  physical_dof_map.dof_indices(elem, physical_dof_indices_var, idx_var);
145  localIdx = physical_param_system.current_solution(physical_dof_indices_var[0]);
146  perf_log.pop("Define physical params");
147 
148  // Restart the FE to the "geometry" of the element
149  // -> Determines quadrature points, shape functions ...
150  // !!! User can change the points to be used (can use other mesh's points
151  // instead of the quadrature points)
152  fe->reinit (elem);
153 
154  perf_log.push("Matrix manipulations");
155  Ke.resize (n_dofs, n_dofs);
156  Fe.resize (n_dofs);
157 
158  // Set the positions of the sub-matrices
159  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
160  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
161  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
162 
163  Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
164  Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
165  Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
166 
167  Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
168  Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
169  Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
170 
171  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
172  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
173  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
174  perf_log.push("Matrix manipulations");
175 
176  perf_log.push("Matrix element calculations");
177  // For each quadrature point determinate the sub-matrices elements
178  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
179  {
180 
181  perf_log.push("Rigidity","Matrix element calculations");
182  // Internal tension
183  alpha_micro = weight_mask.get_alpha_micro(qp_points[qp]);
184 
185  // Internal tension
186  Update_SubK(Kuu, qp, 0, 0, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
187  Update_SubK(Kuv, qp, 0, 1, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
188  Update_SubK(Kuw, qp, 0, 2, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
189 
190  Update_SubK(Kvu, qp, 1, 0, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
191  Update_SubK(Kvv, qp, 1, 1, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
192  Update_SubK(Kvw, qp, 1, 2, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
193 
194  Update_SubK(Kwu, qp, 2, 0, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
195  Update_SubK(Kwv, qp, 2, 1, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
196  Update_SubK(Kww, qp, 2, 2, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha_micro);
197 
198  perf_log.pop("Rigidity","Matrix element calculations");
199  // Gravity
200  // if(z_force)
201  // {
202  for (unsigned int i=0; i<n_w_dofs; i++)
203  {
204  Fw(i) -= alpha_micro * JxW[qp] * phi[i][qp];
205  }
206  }
207  }
208 
209  // Apply constraints
210  perf_log.push("Constraints","Matrix element calculations");
211  dof_map.heterogenously_constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
212  perf_log.pop("Constraints","Matrix element calculations");
213 
214  perf_log.push("Adding elements");
215  system.matrix->add_matrix (Ke, dof_indices);
216  system.rhs->add_vector (Fe, dof_indices);
217  perf_log.pop("Adding elements");
218  }
219 }
220 
221 void assemble_elasticity_anisotropic( libMesh::EquationSystems& es,
222  const std::string& system_name,
223  carl::anisotropic_elasticity_tensor_cubic_sym& anisotropy_obj_input)
224 {
225  libmesh_assert_equal_to (system_name, "Elasticity");
226 
227  libMesh::PerfLog perf_log ("Matrix Assembly (Heterogeneous elasticity)",MASTER_bPerfLog_assemble_fem);
228 
229  perf_log.push("Preamble");
230  // Set up mesh
231  const libMesh::MeshBase& mesh = es.get_mesh();
232 
233  const unsigned int dim = mesh.mesh_dimension();
234 
235  // - Set up physical properties system ------------------------------------
236  int localIdx;
237 
238  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
239  const unsigned int idx_var = physical_param_system.variable_number ("Index");
240 
241  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
242  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
243 
244  // - Set up elasticity system ---------------------------------------------
245  libMesh::LinearImplicitSystem& system = es.get_system<libMesh::LinearImplicitSystem>("Elasticity");
246 
247  const unsigned int n_components = 3;
248  const unsigned int u_var = system.variable_number ("u");
249  const unsigned int v_var = system.variable_number ("v");
250  const unsigned int w_var = system.variable_number ("w");
251 
252  const libMesh::DofMap& dof_map = system.get_dof_map();
253  libMesh::FEType fe_type = dof_map.variable_type(u_var);
254 
255  // Set up pointers to FEBase's of dimension dim and FE type fe_type
256  // -> 3D elements
257  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
258  libMesh::QGauss qrule (dim, fe_type.default_quadrature_order());
259  fe->attach_quadrature_rule (&qrule);
260 
261  // -> Faces
262  libMesh::UniquePtr<libMesh::FEBase> fe_face (libMesh::FEBase::build(dim, fe_type));
263  libMesh::QGauss qface(dim-1, fe_type.default_quadrature_order());
264  fe_face->attach_quadrature_rule (&qface);
265 
266  // Jacobian
267  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
268 
269  // Shape functions
270  const std::vector<std::vector<libMesh::Real> >& phi = fe->get_phi();
271 
272  // Shape functions derivatives
273  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
274 
275  // Quadrature points
276  const std::vector<libMesh::Point>& qp_points = fe->get_xyz();
277 
278  // Weights for the Arlequin method
279  double alpha = 1;
280 
281  libMesh::DenseMatrix<libMesh::Number> Ke;
282  libMesh::DenseVector<libMesh::Number> Fe;
283 
284  libMesh::DenseSubMatrix<libMesh::Number>
285  Kuu(Ke), Kuv(Ke), Kuw(Ke),
286  Kvu(Ke), Kvv(Ke), Kvw(Ke),
287  Kwu(Ke), Kwv(Ke), Kww(Ke);
288 
289  libMesh::DenseSubVector<libMesh::Number>
290  Fu(Fe),
291  Fv(Fe),
292  Fw(Fe);
293 
294  std::vector<libMesh::dof_id_type> dof_indices;
295  std::vector<libMesh::dof_id_type> dof_indices_u;
296  std::vector<libMesh::dof_id_type> dof_indices_v;
297  std::vector<libMesh::dof_id_type> dof_indices_w;
298 
299  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
300  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
301 
302  perf_log.pop("Preamble");
303 
304  // For each element
305  for ( ; el != end_el; ++el)
306  {
307  // Get its pointer
308  const libMesh::Elem* elem = *el;
309 
310  perf_log.push("Define DoF");
311  // The total DoF indices, and those associated to each variable
312  dof_map.dof_indices (elem, dof_indices);
313  dof_map.dof_indices (elem, dof_indices_u, u_var);
314  dof_map.dof_indices (elem, dof_indices_v, v_var);
315  dof_map.dof_indices (elem, dof_indices_w, w_var);
316 
317  const unsigned int n_dofs = dof_indices.size();
318  const unsigned int n_u_dofs = dof_indices_u.size();
319  const unsigned int n_v_dofs = dof_indices_v.size();
320  const unsigned int n_w_dofs = dof_indices_w.size();
321 
322  perf_log.pop("Define DoF");
323 
324  perf_log.push("Define physical params");
325  // The DoF and values of the physical system
326  physical_dof_map.dof_indices(elem, physical_dof_indices_var, idx_var);
327  localIdx = physical_param_system.current_solution(physical_dof_indices_var[0]);
328  perf_log.pop("Define physical params");
329 
330  // Restart the FE to the "geometry" of the element
331  // -> Determines quadrature points, shape functions ...
332  // !!! User can change the points to be used (can use other mesh's points
333  // instead of the quadrature points)
334  fe->reinit (elem);
335 
336  perf_log.push("Matrix manipulations");
337  Ke.resize (n_dofs, n_dofs);
338  Fe.resize (n_dofs);
339 
340  // Set the positions of the sub-matrices
341  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
342  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
343  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
344 
345  Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
346  Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
347  Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
348 
349  Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
350  Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
351  Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
352 
353  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
354  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
355  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
356  perf_log.push("Matrix manipulations");
357 
358  perf_log.push("Matrix element calculations");
359  // For each quadrature point determinate the sub-matrices elements
360  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
361  {
362 
363  perf_log.push("Rigidity","Matrix element calculations");
364  // Internal tension
365  alpha = 1;
366 
367  // Internal tension
368  Update_SubK(Kuu, qp, 0, 0, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
369  Update_SubK(Kuv, qp, 0, 1, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
370  Update_SubK(Kuw, qp, 0, 2, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
371 
372  Update_SubK(Kvu, qp, 1, 0, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
373  Update_SubK(Kvv, qp, 1, 1, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
374  Update_SubK(Kvw, qp, 1, 2, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
375 
376  Update_SubK(Kwu, qp, 2, 0, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
377  Update_SubK(Kwv, qp, 2, 1, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
378  Update_SubK(Kww, qp, 2, 2, dphi, n_components, n_u_dofs, JxW, anisotropy_obj_input, localIdx, alpha);
379 
380  perf_log.pop("Rigidity","Matrix element calculations");
381  // Gravity
382  // if(z_force)
383  // {
384  for (unsigned int i=0; i<n_w_dofs; i++)
385  {
386  Fw(i) -= alpha * JxW[qp] * phi[i][qp];
387  }
388  }
389  }
390 
391  // Apply constraints
392  perf_log.push("Constraints","Matrix element calculations");
393  dof_map.heterogenously_constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
394  perf_log.pop("Constraints","Matrix element calculations");
395 
396  perf_log.push("Adding elements");
397  system.matrix->add_matrix (Ke, dof_indices);
398  system.rhs->add_vector (Fe, dof_indices);
399  perf_log.pop("Adding elements");
400  }
401 }
Update_SubK
void Update_SubK(libMesh::DenseSubMatrix< libMesh::Number > &SubK, unsigned int qp, unsigned int C_i, unsigned int C_k, const std::vector< std::vector< libMesh::RealGradient > > &dphi, const unsigned int n_components, const unsigned int n_u_dofs, const std::vector< libMesh::Real > &JxW, carl::anisotropic_elasticity_tensor_cubic_sym &anisotropy_obj_input, int grain_idx, double cte)
assemble_functions_elasticity_anisotropy_3D.h
assemble_elasticity_anisotropic_with_weight
void assemble_elasticity_anisotropic_with_weight(libMesh::EquationSystems &es, const std::string &system_name, carl::weight_parameter_function &weight_mask, carl::anisotropic_elasticity_tensor_cubic_sym &anisotropy_obj_input)
Definition: assemble_functions_elasticity_anisotropy_3D.cpp:38
carl::weight_parameter_function::get_alpha_micro
double get_alpha_micro(const libMesh::Point &qpoint)
Definition: weight_parameter_function.h:100
assemble_elasticity_anisotropic
void assemble_elasticity_anisotropic(libMesh::EquationSystems &es, const std::string &system_name, carl::anisotropic_elasticity_tensor_cubic_sym &anisotropy_obj_input)
Definition: assemble_functions_elasticity_anisotropy_3D.cpp:221
carl::weight_parameter_function
Definition: weight_parameter_function.h:16
Update_SubK_anisotropic
void Update_SubK_anisotropic(libMesh::DenseSubMatrix< libMesh::Number > &SubK, unsigned int qp, unsigned int C_i, unsigned int C_k, const std::vector< std::vector< libMesh::RealGradient > > &dphi, const unsigned int n_components, const unsigned int n_u_dofs, const std::vector< libMesh::Real > &JxW, carl::anisotropic_elasticity_tensor_cubic_sym &anisotropy_obj_input, int grain_idx, double cte)
Definition: assemble_functions_elasticity_anisotropy_3D.cpp:10
MASTER_bPerfLog_assemble_fem
const bool MASTER_bPerfLog_assemble_fem
Definition: assemble_functions_elasticity_3D.h:17