MAST
level_set_eigenproblem_assembly.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 // MAST includes
23 #include "level_set/sub_cell_fe.h"
26 #include "base/nonlinear_system.h"
29 #include "base/elem_base.h"
30 #include "numerics/utility.h"
31 #include "mesh/geom_elem.h"
32 
33 // libMesh includes
34 #include "libmesh/nonlinear_solver.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/sparse_matrix.h"
37 #include "libmesh/dof_map.h"
38 
39 
40 
44 _level_set (nullptr),
45 _intersection (nullptr),
46 _velocity (nullptr) {
47 
48 }
49 
50 
51 
53 
54  if (_intersection)
55  delete _intersection;
56 }
57 
58 
61 
62  libmesh_assert(_level_set);
63  return *_intersection;
64 }
65 
66 
67 void
70 
71  libmesh_assert(!_level_set);
72  libmesh_assert(!_intersection);
73  libmesh_assert(_system);
74 
75  _level_set = &level_set;
77 }
78 
79 
80 
81 void
83 
84  _level_set = nullptr;
85 
86  if (_intersection) {
87  delete _intersection;
88  _intersection = nullptr;
89  }
90 }
91 
92 
93 
94 
95 void
98 
99  libmesh_assert(_level_set);
100  libmesh_assert(_intersection);
101  libmesh_assert(_system);
102  libmesh_assert(!_velocity);
103 
104  _velocity = &velocity;
105 }
106 
107 
108 
109 void
111 
112  _velocity = nullptr;
113 }
114 
115 
116 
117 void
119 eigenproblem_assemble(libMesh::SparseMatrix<Real>* A,
120  libMesh::SparseMatrix<Real>* B) {
121 
122  libmesh_assert(_system);
123  libmesh_assert(_discipline);
124  libmesh_assert(_elem_ops);
125  libmesh_assert(_level_set);
126 
127  MAST::NonlinearSystem& eigen_sys =
128  dynamic_cast<MAST::NonlinearSystem&>(_system->system());
129 
130  libMesh::SparseMatrix<Real>
131  &matrix_A = *A,
132  &matrix_B = *B;
133 
134  matrix_A.zero();
135  matrix_B.zero();
136 
137  // build localized solutions if needed
138  std::unique_ptr<libMesh::NumericVector<Real> >
139  localized_solution;
140 
141  if (_base_sol)
142  localized_solution.reset(build_localized_vector(eigen_sys,
143  *_base_sol).release());
144 
145 
146  // iterate over each element, initialize it and get the relevant
147  // analysis quantities
148  RealVectorX sol;
149  RealMatrixX mat_A, mat_B;
150  std::vector<libMesh::dof_id_type> dof_indices;
151  const libMesh::DofMap& dof_map = eigen_sys.get_dof_map();
152 
153 
154  libMesh::MeshBase::const_element_iterator el =
155  eigen_sys.get_mesh().active_local_elements_begin();
156  const libMesh::MeshBase::const_element_iterator end_el =
157  eigen_sys.get_mesh().active_local_elements_end();
158 
160  &ops = dynamic_cast<MAST::EigenproblemAssemblyElemOperations&>(*_elem_ops);
161 
162  for ( ; el != end_el; ++el) {
163 
164  const libMesh::Elem* elem = *el;
165 
166  _intersection->init(*_level_set, *elem, eigen_sys.time,
167  eigen_sys.get_mesh().max_elem_id(),
168  eigen_sys.get_mesh().max_node_id());
169 
170  dof_map.dof_indices (elem, dof_indices);
171 
173 
174  // get the solution
175  unsigned int ndofs = (unsigned int)dof_indices.size();
176 
177  sol.setZero(ndofs);
178  mat_A.setZero(ndofs, ndofs);
179  mat_B.setZero(ndofs, ndofs);
180 
181  // if the base solution is provided, then tell the element about it
182  if (_base_sol) {
183  for (unsigned int i=0; i<dof_indices.size(); i++)
184  sol(i) = (*localized_solution)(dof_indices[i]);
185  }
186 
187  /*const std::vector<const libMesh::Elem *> &
188  elems_hi = _intersection->get_sub_elems_positive_phi();
189 
190 
191  std::vector<const libMesh::Elem*>::const_iterator
192  hi_sub_elem_it = elems_hi.begin(),
193  hi_sub_elem_end = elems_hi.end();
194 
195  for (; hi_sub_elem_it != hi_sub_elem_end; hi_sub_elem_it++ )*/ {
196 
197  //const libMesh::Elem* sub_elem = *hi_sub_elem_it;
198 
199  MAST::GeomElem geom_elem;
200  ops.set_elem_data(elem->dim(), *elem, geom_elem);
201  geom_elem.init(*elem, *_system);
202 
203  ops.init(geom_elem);
204  ops.set_elem_solution(sol);
205  ops.elem_calculations(mat_A, mat_B);
206  ops.clear_elem();
207 
208  // copy to the libMesh matrix for further processing
209  DenseRealMatrix A, B;
210  MAST::copy(A, mat_A);
211  MAST::copy(B, mat_B);
212 
213  // constrain the element matrices.
214  dof_map.constrain_element_matrix(A, dof_indices);
215  dof_map.constrain_element_matrix(B, dof_indices);
216 
217  matrix_A.add_matrix (A, dof_indices); // load independent
218  matrix_B.add_matrix (B, dof_indices); // load dependent
219  }
220  }
221 
222  _intersection->clear();
223  }
224 
225  // finalize the data structures
226  A->close();
227  B->close();
228 }
229 
230 
231 
232 
233 bool
236  libMesh::SparseMatrix<Real>* sensitivity_A,
237  libMesh::SparseMatrix<Real>* sensitivity_B) {
238 
239  libmesh_assert(_system);
240  libmesh_assert(_discipline);
241  libmesh_assert(_elem_ops);
242  // we need the velocity for topology parameter
243  if (f.is_topology_parameter()) libmesh_assert(_velocity);
244 
245  MAST::NonlinearSystem& eigen_sys =
246  dynamic_cast<MAST::NonlinearSystem&>(_system->system());
247 
248  libMesh::SparseMatrix<Real>& matrix_A = *sensitivity_A;
249  libMesh::SparseMatrix<Real>& matrix_B = *sensitivity_B;
250 
251  matrix_A.zero();
252 
253  // build localized solutions if needed
254  std::unique_ptr<libMesh::NumericVector<Real> >
255  localized_solution,
256  localized_solution_sens;
257 
258  if (_base_sol) {
259 
260  localized_solution.reset(build_localized_vector(eigen_sys,
261  *_base_sol).release());
262 
263  // make sure that the sensitivity was also provided
264  libmesh_assert(_base_sol_sensitivity);
265  localized_solution_sens.reset(build_localized_vector(eigen_sys,
266  *_base_sol_sensitivity).release());
267  }
268 
269 
270  // iterate over each element, initialize it and get the relevant
271  // analysis quantities
272  RealVectorX sol, dsol;
273  RealMatrixX mat_A, mat_B, mat2_A, mat2_B;
274  std::vector<libMesh::dof_id_type> dof_indices;
275  const libMesh::DofMap& dof_map = eigen_sys.get_dof_map();
276 
277 
278  libMesh::MeshBase::const_element_iterator el =
279  eigen_sys.get_mesh().active_local_elements_begin();
280  const libMesh::MeshBase::const_element_iterator end_el =
281  eigen_sys.get_mesh().active_local_elements_end();
282 
284  &ops = dynamic_cast<MAST::EigenproblemAssemblyElemOperations&>(*_elem_ops);
285 
286  for ( ; el != end_el; ++el) {
287 
288  const libMesh::Elem* elem = *el;
289 
290  _intersection->init(*_level_set, *elem, eigen_sys.time,
291  eigen_sys.get_mesh().max_elem_id(),
292  eigen_sys.get_mesh().max_node_id());
293 
294 
296 
297  dof_map.dof_indices (elem, dof_indices);
298 
299  // get the solution
300  unsigned int ndofs = (unsigned int)dof_indices.size();
301  dsol.setZero(ndofs);
302  sol.setZero(ndofs);
303  mat_A.setZero(ndofs, ndofs);
304  mat_B.setZero(ndofs, ndofs);
305  mat2_A.setZero(ndofs, ndofs);
306  mat2_B.setZero(ndofs, ndofs);
307 
308  // if the base solution is provided, tell the element about it
309  if (_base_sol) {
310 
311  // set the element's base solution
312  for (unsigned int i=0; i<dof_indices.size(); i++) {
313 
314  sol(i) = (*localized_solution)(dof_indices[i]);
315  dsol(i) = (*localized_solution_sens)(dof_indices[i]);
316  }
317  }
318 
319  const std::vector<const libMesh::Elem *> &
321 
322  std::vector<const libMesh::Elem*>::const_iterator
323  hi_sub_elem_it = elems_hi.begin(),
324  hi_sub_elem_end = elems_hi.end();
325 
326  for (; hi_sub_elem_it != hi_sub_elem_end; hi_sub_elem_it++ ) {
327 
328  const libMesh::Elem* sub_elem = *hi_sub_elem_it;
329 
331  ops.set_elem_data(elem->dim(), *elem, geom_elem);
332  geom_elem.init(*sub_elem, *_system, *_intersection);
333 
334  ops.init(geom_elem);
335  ops.set_elem_solution(sol);
336  if (_base_sol)
338 
340  _base_sol!=nullptr,
341  mat_A,
342  mat_B);
343 
344  if (f.is_topology_parameter()) {
346  _base_sol!=nullptr,
347  *_velocity,
348  mat2_A,
349  mat2_B);
350 
351  mat_A += mat2_A;
352  mat_B += mat2_B;
353  }
354  ops.clear_elem();
355 
356  // copy to the libMesh matrix for further processing
357  DenseRealMatrix A, B;
358  MAST::copy(A, mat_A);
359  MAST::copy(B, mat_B);
360 
361  // constrain the element matrices.
362  dof_map.constrain_element_matrix(A, dof_indices);
363  dof_map.constrain_element_matrix(B, dof_indices);
364 
365  matrix_A.add_matrix (A, dof_indices);
366  matrix_B.add_matrix (B, dof_indices);
367  }
368  }
369  _intersection->clear();
370  }
371 
372  // finalize the data structures
373  sensitivity_A->close();
374  sensitivity_B->close();
375 
376  return true;
377 }
378 
379 
MAST::AssemblyElemOperations * _elem_ops
provides assembly elem operations for use by this class
MAST::NonlinearSystem & system()
libMesh::DenseMatrix< Real > DenseRealMatrix
This class implements a system for solution of nonlinear systems.
virtual void elem_topology_sensitivity_calculations(const MAST::FunctionBase &f, bool base_sol, const MAST::FieldFunction< RealVectorX > &vel, RealMatrixX &mat_A, RealMatrixX &mat_B)=0
performs the element topology sensitivity calculations over elem.
virtual void set_elem_solution_sensitivity(const RealVectorX &sol)
sets the element solution sensitivity
MAST::FieldFunction< RealVectorX > * _velocity
virtual void init(const libMesh::Elem &elem, const MAST::SystemInitialization &sys_init)
This method should not get called for this class.
virtual void elem_calculations(RealMatrixX &mat_A, RealMatrixX &mat_B)=0
performs the element calculations over elem, and returns the element matrices for the eigenproblem ...
virtual void set_elem_solution(const RealVectorX &sol)
sets the element solution
virtual void set_level_set_function(MAST::FieldFunction< Real > &level_set)
attaches level set function to this
void init(const MAST::FieldFunction< Real > &phi, const libMesh::Elem &e, const Real t, unsigned int max_elem_id, unsigned int max_node_id)
MAST::SystemInitialization * _system
System for which this assembly is performed.
MAST::PhysicsDisciplineBase * _discipline
PhysicsDisciplineBase object for which this class is assembling.
virtual void set_elem_data(unsigned int dim, const libMesh::Elem &ref_elem, MAST::GeomElem &elem) const =0
some analyses may want to set additional element data before initialization of the GeomElem...
virtual void clear_level_set_function()
clears association with level set function
Assembles the system of equations for an eigenproblem of type .
Matrix< Real, Dynamic, Dynamic > RealMatrixX
const libMesh::NumericVector< Real > * _base_sol
base solution about which this eigenproblem is defined.
void copy(DenseRealMatrix &m1, const RealMatrixX &m2)
Definition: utility.h:167
This class inherits from MAST::GeomElem and provides an interface to initialize FE objects on sub-ele...
Matrix< Real, Dynamic, 1 > RealVectorX
virtual void init(const MAST::GeomElem &elem)=0
initializes the object for calculation of element quantities for the specified elem.
virtual ~LevelSetEigenproblemAssembly()
destructor resets the association of this assembly object with the system
virtual void set_level_set_velocity_function(MAST::FieldFunction< RealVectorX > &velocity)
the velocity function used to calculate topology sensitivity
std::unique_ptr< libMesh::NumericVector< Real > > build_localized_vector(const libMesh::System &sys, const libMesh::NumericVector< Real > &global) const
localizes the parallel vector so that the local copy stores all values necessary for calculation of t...
virtual void eigenproblem_assemble(libMesh::SparseMatrix< Real > *A, libMesh::SparseMatrix< Real > *B)
assembles the matrices for eigenproblem depending on the analysis type
This class acts as a wrapper around libMesh::Elem for the purpose of providing a uniform interface fo...
Definition: geom_elem.h:59
const libMesh::NumericVector< Real > * _base_sol_sensitivity
sensitivity of base solution may be needed for sensitivity analysis.
const std::vector< const libMesh::Elem * > & get_sub_elems_positive_phi() const
virtual void elem_sensitivity_calculations(const MAST::FunctionBase &f, bool base_sol, RealMatrixX &mat_A, RealMatrixX &mat_B)=0
performs the element sensitivity calculations over elem, and returns the element matrices for the eig...
virtual void clear_elem()
clears the element initialization
virtual void init(const libMesh::Elem &elem, const MAST::SystemInitialization &sys_init)
initialize the object for the specified reference elem.
Definition: geom_elem.cpp:128
virtual bool eigenproblem_sensitivity_assemble(const MAST::FunctionBase &f, libMesh::SparseMatrix< Real > *sensitivity_A, libMesh::SparseMatrix< Real > *sensitivity_B)
Assembly function.
LevelSetEigenproblemAssembly()
constructor associates this assembly object with the system
virtual bool is_topology_parameter() const
Definition: function_base.h:97
void clear()
clears the data structures
virtual void clear_level_set_velocity_function()
clears the velocity function