MAST
interface_dof_handler.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 
21 // MAST includes
27 #include "base/nonlinear_system.h"
28 
29 // libMesh includes
30 #include "libmesh/dof_map.h"
31 
32 
34 _sys_init (nullptr),
35 _phi (nullptr) {
36 
37 }
38 
39 
40 
42 
43 }
44 
45 
46 void
48  MAST::LevelSetIntersection& intersection,
50 
51  libmesh_assert(!_sys_init);
52 
53  _sys_init = &sys_init;
54  _phi = φ
55 
56  const MAST::NonlinearSystem &system = sys_init.system();
57  const libMesh::DofMap &dof_map = system.get_dof_map();
58  const libMesh::MeshBase &mesh = system.get_mesh();
59 
60  std::set<const libMesh::Node*> negative_phi_nodes;
61  std::vector<const libMesh::Elem*>
62  elems_to_factor,
63  negative_phi_elems;
64  elems_to_factor.reserve(mesh.n_local_elem()),
65  negative_phi_elems.reserve(mesh.n_local_elem());
66 
67  // iterate over all the elements identify if elements on the interface
68  // need to factor their stiffness matrices
69  libMesh::MeshBase::const_element_iterator
70  el = mesh.active_local_elements_begin(),
71  end_el = mesh.active_local_elements_end();
72 
73 
74  for ( ; el != end_el; ++el) {
75 
76  const libMesh::Elem* elem = *el;
77  intersection.init(phi, *elem, system.time,
78  system.get_mesh().max_elem_id(),
79  system.get_mesh().max_node_id());
80  if (intersection.if_elem_has_boundary()) {
81 
82  intersection.get_corner_nodes_on_negative_phi(negative_phi_nodes);
83 
84  std::set<const libMesh::Node*>::const_iterator
85  nd_it = negative_phi_nodes.begin(),
86  nd_end = negative_phi_nodes.end();
87 
88  for ( ; nd_it != nd_end; nd_it++) {
89 
90  MAST::MaterialPatch patch;
91  patch.init(*elem, **nd_it, phi, system.time);
92  std::set<const libMesh::Elem*>
93  patch_elems_to_factor = patch.get_elems_to_factor();
94  std::set<const libMesh::Elem*>::const_iterator
95  el_it = patch_elems_to_factor.begin(),
96  el_end = patch_elems_to_factor.end();
97 
98  for ( ; el_it != el_end; el_it++)
99  elems_to_factor.push_back(*el_it);
100  }
101  }
102 
103  negative_phi_nodes.clear();
104  intersection.clear();
105  }
106 
107  // now create a unique set from the vector of elements to facotr
108  std::set<const libMesh::Elem*> elems_to_factor_set(elems_to_factor.begin(),
109  elems_to_factor.end());
110 
111  elems_to_factor.clear();
112 
113  // now, create an entry in the map for each one of these elements
114  std::set<const libMesh::Elem*>::const_iterator
115  el_it = elems_to_factor_set.begin(),
116  el_end = elems_to_factor_set.end();
117 
118  std::vector<libMesh::dof_id_type> dof_ids;
119  std::pair<const libMesh::Elem*, RealVectorX> elem_sol_pair;
120 
121  for ( ; el_it != el_end; el_it++) {
122 
123  dof_map.dof_indices(*el_it, dof_ids);
124  elem_sol_pair.first = *el_it;
125  elem_sol_pair.second = RealVectorX::Zero(dof_ids.size());
126  _elem_sol.insert(elem_sol_pair);
127  }
128 }
129 
130 
131 bool
133 
134  libmesh_assert(_sys_init);
135 
136  return _elem_sol.count(&elem);
137 }
138 
139 
140 void
142 partition_local_elem_rows(const libMesh::Elem& elem,
143  std::vector<libMesh::dof_id_type>& material_dofs,
144  std::vector<libMesh::dof_id_type>& void_dofs) {
145 
146 
147  const MAST::NonlinearSystem &system = _sys_init->system();
148  Real val = 0.;
149  unsigned int
150  nvars = _sys_init->n_vars(),
151  n_nodes = elem.n_nodes();
152 
153  // currently only implemented for quad4 and Lagrange functions
154  libmesh_assert_equal_to(elem.type(), libMesh::QUAD4);
155  // the number of dofs is assumed to be a multiple of the number of variables.
156  // currently, it is assumed that the Lagrange shape functions are used,
157  // so that each node has dofs for each variable.
158  libmesh_assert_equal_to(system.variable_type(0).family, libMesh::LAGRANGE);
159 
160  // the dofs are assumed to be sequenced in variable major format
161  // So, we loop over the nodes and if a node is in the void domain, all
162  // dofs on that will be factored out.
163  material_dofs.reserve(nvars*elem.n_nodes());
164  void_dofs.reserve(nvars*elem.n_nodes());
165 
166  for (unsigned int i=0; i<elem.n_nodes(); i++) {
167 
168  (*_phi)(*elem.node_ptr(i), system.time, val);
169  if (val < 0.) {
170 
171  // mark all dofs for this node to be in void
172  for (unsigned int j=0; j<nvars; j++)
173  void_dofs.push_back(j*n_nodes+i);
174  }
175  else {
176 
177  // mark all dofs for this node to be in material
178  for (unsigned int j=0; j<nvars; j++)
179  material_dofs.push_back(j*n_nodes+i);
180  }
181  }
182 
183  // now sort the dofs
184  std::sort(void_dofs.begin(), void_dofs.end());
185  std::sort(material_dofs.begin(), material_dofs.end());
186 }
187 
188 
189 void
191 partition_global_elem_rows(const libMesh::Elem& elem,
192  std::vector<libMesh::dof_id_type>& material_dofs,
193  std::vector<libMesh::dof_id_type>& void_dofs) {
194 
195  // get the local dof partitioning
196  this->partition_local_elem_rows(elem, material_dofs, void_dofs);
197 
198  //
199  // now replace the dof ids with the global ids
200  //
201  const MAST::NonlinearSystem &system = _sys_init->system();
202  const libMesh::DofMap &dof_map = system.get_dof_map();
203  std::vector<libMesh::dof_id_type>
204  dof_ids;
205 
206  dof_map.dof_indices(&elem, dof_ids);
207 
208  for (unsigned int i=0; i<material_dofs.size(); i++)
209  material_dofs[i] = dof_ids[material_dofs[i]];
210 
211  for (unsigned int i=0; i<void_dofs.size(); i++)
212  void_dofs[i] = dof_ids[void_dofs[i]];
213 }
214 
215 
216 void
218  RealVectorX& elem_sol) {
219 
220  // make sure that the element has been identified for factorization
221  std::map<const libMesh::Elem*, RealVectorX>::const_iterator
222  it = _elem_sol.find(&elem),
223  end = _elem_sol.end();
224 
225  libmesh_assert(it != end);
226 
227  std::vector<libMesh::dof_id_type>
228  material_dofs,
229  void_dofs;
230 
231  const RealVectorX
232  &free_sol = it->second;
233 
234  this->partition_local_elem_rows(elem, material_dofs, void_dofs);
235 
236  // overwirte the void dofs from the stored data
237  for (unsigned int i=0; i<void_dofs.size(); i++)
238  elem_sol(void_dofs[i]) = free_sol(void_dofs[i]);
239 }
240 
241 
242 
243 void
245 element_factored_jacobian(const libMesh::Elem& elem,
246  const RealMatrixX& jac,
247  std::vector<libMesh::dof_id_type>& material_dof_ids,
248  RealMatrixX& jac_factored_uu) {
249 
250  std::vector<libMesh::dof_id_type>
251  void_dof_ids;
252 
254  jac_uu,
255  jac_uf,
256  jac_fu,
257  jac_ff;
258 
259  this->element_factored_jacobian(elem,
260  jac,
261  material_dof_ids,
262  void_dof_ids,
263  jac_uu,
264  jac_uf,
265  jac_fu,
266  jac_ff,
267  jac_factored_uu);
268 
269 }
270 
271 
272 
273 void
275 element_factored_residual_and_jacobian(const libMesh::Elem& elem,
276  const RealMatrixX& jac,
277  const RealVectorX& res,
278  std::vector<libMesh::dof_id_type>& material_dof_ids,
279  RealMatrixX& jac_factored_uu,
280  RealVectorX& res_factored_u) {
281 
282  std::vector<libMesh::dof_id_type>
283  void_dof_ids;
284 
286  jac_uu,
287  jac_uf,
288  jac_fu,
289  jac_ff,
290  jac_ff_inv;
291 
292  this->element_factored_jacobian(elem,
293  jac,
294  material_dof_ids,
295  void_dof_ids,
296  jac_uu,
297  jac_uf,
298  jac_fu,
299  jac_ff,
300  jac_factored_uu);
301 
302  unsigned int
303  nu = material_dof_ids.size(),
304  nf = void_dof_ids.size();
305 
307  res_f = RealVectorX::Zero(nf);
308 
309  res_factored_u = RealVectorX::Zero(nu);
310 
311  for (unsigned int i=0; i<material_dof_ids.size(); i++)
312  res_factored_u(i) = res(material_dof_ids[i]);
313 
314  for (unsigned int i=0; i<void_dof_ids.size(); i++)
315  res_f(i) = res(void_dof_ids[i]);
316 
317  _compute_matrix_inverse(jac_ff, jac_ff_inv);
318 
319  res_factored_u -= jac_uf * jac_ff_inv * res_f;
320 }
321 
322 
323 
324 void
326 element_factored_jacobian(const libMesh::Elem& elem,
327  const RealMatrixX& jac,
328  std::vector<libMesh::dof_id_type>& material_dof_ids,
329  std::vector<libMesh::dof_id_type>& void_dof_ids,
330  RealMatrixX& jac_uu,
331  RealMatrixX& jac_uf,
332  RealMatrixX& jac_fu,
333  RealMatrixX& jac_ff,
334  RealMatrixX& jac_factored_uu) {
335 
336 
337  material_dof_ids.clear();
338  void_dof_ids.clear();
339 
340  this->partition_local_elem_rows(elem, material_dof_ids, void_dof_ids);
341 
342  unsigned int
343  nu = material_dof_ids.size(),
344  nf = void_dof_ids.size();
345 
346  jac_uu = RealMatrixX::Zero(nu, nu);
347  jac_uf = RealMatrixX::Zero(nu, nf);
348  jac_fu = RealMatrixX::Zero(nf, nu);
349  jac_ff = RealMatrixX::Zero(nf, nf);
350 
352  jac_ff_inv;
353 
354  //
355  // partition the matrices
356  //
357  for (unsigned int i=0; i<nu; i++) {
358  for (unsigned int j=0; j<nu; j++)
359  jac_uu(i,j) = jac(material_dof_ids[i], material_dof_ids[j]);
360 
361  for (unsigned int j=0; j<nf; j++)
362  jac_uf(i,j) = jac(material_dof_ids[i], void_dof_ids[j]);
363  }
364 
365  for (unsigned int i=0; i<nf; i++) {
366  for (unsigned int j=0; j<nu; j++)
367  jac_fu(i,j) = jac(void_dof_ids[i], material_dof_ids[j]);
368 
369  for (unsigned int j=0; j<nf; j++)
370  jac_ff(i,j) = jac(void_dof_ids[i], void_dof_ids[j]);
371  }
372 
373  _compute_matrix_inverse(jac_ff, jac_ff_inv);
374 
375  jac_factored_uu = jac_uu - jac_uf * jac_ff_inv * jac_fu;
376 }
377 
378 
379 void
381 update_factored_element_solution(const libMesh::Elem& elem,
382  const RealMatrixX& res,
383  const RealMatrixX& jac,
384  const RealMatrixX& sol,
385  const RealMatrixX& dsol,
386  RealVectorX& updated_sol) {
387 
388  std::map<const libMesh::Elem*, RealVectorX>::iterator
389  it = _elem_sol.find(&elem),
390  end = _elem_sol.end();
391 
392  libmesh_assert(it != end);
393 
394  //
395  // [Juu Juf] {dxu} = - {f_u}
396  // [Jfu Jff] {dxf} = - {f_f}
397  // so,
398  // {dxf} = inv(Jff) (-{f_f} - Jfu {dxu})
399  //
400 
401  std::vector<libMesh::dof_id_type>
402  material_dof_ids,
403  void_dof_ids;
404 
405  this->partition_local_elem_rows(elem, material_dof_ids, void_dof_ids);
406 
407  // we should not be capturing elements that are purely in void or material
408  libmesh_assert_greater(material_dof_ids.size(), 0);
409  libmesh_assert_greater( void_dof_ids.size(), 0);
410 
411  unsigned int
412  nu = material_dof_ids.size(),
413  nf = void_dof_ids.size();
414 
416  jac_fu = RealMatrixX::Zero(nf, nu),
417  jac_ff = RealMatrixX::Zero(nf, nf),
418  jac_ff_inv;
419 
421  f_u = RealVectorX::Zero(nu),
422  f_f = RealVectorX::Zero(nf),
423  dx_u = RealVectorX::Zero(nu),
424  dx_f = RealVectorX::Zero(nf);
425 
426 
427  //
428  // partition the matrices
429  //
430  for (unsigned int i=0; i<nf; i++) {
431 
432  for (unsigned int j=0; j<nu; j++)
433  jac_fu(i,j) = jac(void_dof_ids[i], material_dof_ids[j]);
434 
435  for (unsigned int j=0; j<nf; j++)
436  jac_ff(i,j) = jac(void_dof_ids[i], void_dof_ids[j]);
437 
438  f_f(i) = res(void_dof_ids[i]);
439  }
440 
441  for (unsigned int i=0; i<nu; i++)
442  dx_u(i)= dsol(material_dof_ids[i]);
443 
444  _compute_matrix_inverse(jac_ff, jac_ff_inv);
445 
446  dx_f = jac_ff_inv * (-f_f - jac_fu * dx_u);
447 
448  updated_sol = sol;
449  for (unsigned int i=0; i<nf; i++)
450  updated_sol(void_dof_ids[i]) += dx_f(i);
451 
452 
453  //
454  // now update the solution in the map
455  //
456  it->second = updated_sol;
457 }
458 
459 
460 void
462  RealMatrixX& mat_inv) {
463 
464  const unsigned int
465  m = mat.rows(),
466  n = mat.cols();
467 
468  libmesh_assert_equal_to(m, n);
469 
470  Eigen::FullPivLU<RealMatrixX> solver(mat);
472  rhs = RealVectorX::Zero(m);
473 
474  mat_inv = RealMatrixX::Zero(m, m);
475 
476  for (unsigned int i=0; i<m; i++) {
477 
478  rhs = RealVectorX::Zero(m);
479  rhs(i) = 1.;
480 
481  mat_inv.col(i) = solver.solve(rhs);
482  }
483 }
484 
void element_factored_jacobian(const libMesh::Elem &elem, const RealMatrixX &jac, std::vector< libMesh::dof_id_type > &material_dof_ids, RealMatrixX &jac_factored_uu)
a wrapper around the second element_factored_jacobian.
void update_factored_element_solution(const libMesh::Elem &elem, const RealMatrixX &res, const RealMatrixX &jac, const RealMatrixX &sol, const RealMatrixX &dsol, RealVectorX &updated_sol)
MAST::NonlinearSystem & system()
std::map< const libMesh::Elem *, RealVectorX > _elem_sol
This class implements a system for solution of nonlinear systems.
void solution_of_factored_element(const libMesh::Elem &elem, RealVectorX &elem_sol)
updates the components of the solution vector in elem_sol for the void domain using the stored soluti...
void init(const libMesh::Elem &elem, const libMesh::Node &node, const MAST::FieldFunction< Real > &phi, const Real t)
initialize the patch around node of elem.
const std::set< const libMesh::Elem * > & get_elems_to_factor() const
libMesh::Real Real
void init(const MAST::FieldFunction< Real > &phi, const libMesh::Elem &e, const Real t, unsigned int max_elem_id, unsigned int max_node_id)
void _compute_matrix_inverse(const RealMatrixX &mat, RealMatrixX &mat_inv)
MAST::FieldFunction< Real > * _phi
Matrix< Real, Dynamic, Dynamic > RealMatrixX
A patch is defines as the set of elements sharing a node.
void element_factored_residual_and_jacobian(const libMesh::Elem &elem, const RealMatrixX &jac, const RealVectorX &res, std::vector< libMesh::dof_id_type > &material_dof_ids, RealMatrixX &jac_factored_uu, RealVectorX &res_factored_u)
factorizes the residual and jacobian into the components for the dofs on material nodes...
Matrix< Real, Dynamic, 1 > RealVectorX
void get_corner_nodes_on_negative_phi(std::set< const libMesh::Node * > &nodes) const
void partition_global_elem_rows(const libMesh::Elem &elem, std::vector< libMesh::dof_id_type > &material_dofs, std::vector< libMesh::dof_id_type > &void_dofs)
fills the material_dofs and void_dofs with the dofs_ids in the global system corresponding to these d...
void init(const MAST::SystemInitialization &sys_init, MAST::LevelSetIntersection &intersection, MAST::FieldFunction< Real > &phi)
const MAST::SystemInitialization * _sys_init
bool if_factor_element(const libMesh::Elem &elem) const
void partition_local_elem_rows(const libMesh::Elem &elem, std::vector< libMesh::dof_id_type > &material_dofs, std::vector< libMesh::dof_id_type > &void_dofs)
identifies which rows in the element residual vector and rows/columns in the jacobian matrix correspo...
void clear()
clears the data structures