MAST
structural_fluid_interaction_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 
21 // MAST includes
27 #include "base/nonlinear_system.h"
29 #include "numerics/utility.h"
30 #include "mesh/geom_elem.h"
31 
32 // libMesh includes
33 #include "libmesh/numeric_vector.h"
34 #include "libmesh/sparse_matrix.h"
35 #include "libmesh/dof_map.h"
36 #include "libmesh/petsc_nonlinear_solver.h"
37 #include "libmesh/petsc_vector.h"
38 
39 
40 
44 _base_sol (nullptr),
45 _base_sol_sensitivity (nullptr) {
46 
47 }
48 
49 
50 
53 
54 }
55 
56 
57 
58 
59 
60 void
63 
64  _base_sol = nullptr;
65  _base_sol_sensitivity = nullptr;
66 
68 }
69 
70 
71 
72 
73 void
74 MAST::StructuralFluidInteractionAssembly::set_base_solution(const libMesh::NumericVector<Real>& sol,
75  bool if_sens) {
76 
77  if (!if_sens) {
78 
79  // make sure that the pointer has been cleared
80  libmesh_assert(!_base_sol);
81  _base_sol = &sol;
82  }
83  else {
84 
85  libmesh_assert(!_base_sol_sensitivity);
86  _base_sol_sensitivity = &sol;
87  }
88 }
89 
90 
91 
92 
93 void
95 
96  if (!if_sens)
97  _base_sol = nullptr;
98  else
99  _base_sol_sensitivity = nullptr;
100 }
101 
102 
103 
104 
105 void
108 (std::vector<libMesh::NumericVector<Real>*>& basis,
109  std::map<MAST::StructuralQuantityType, RealMatrixX*>& mat_qty_map) {
110 
111  libmesh_assert(_elem_ops);
112 
113  MAST::NonlinearSystem& nonlin_sys = _system->system();
114 
115  unsigned int
116  n_basis = (unsigned int)basis.size();
117 
118  // initialize the quantities to zero matrices
119  std::map<MAST::StructuralQuantityType, RealMatrixX*>::iterator
120  it = mat_qty_map.begin(),
121  end = mat_qty_map.end();
122 
123  for ( ; it != end; it++)
124  *it->second = RealMatrixX::Zero(n_basis, n_basis);
125 
126  // iterate over each element, initialize it and get the relevant
127  // analysis quantities
128  RealVectorX vec, sol;
129  RealMatrixX mat, basis_mat;
130 
131  std::vector<libMesh::dof_id_type> dof_indices;
132  const libMesh::DofMap& dof_map = nonlin_sys.get_dof_map();
133 
134 
135  std::unique_ptr<libMesh::NumericVector<Real> > localized_solution;
136  if (_base_sol)
137  localized_solution.reset(build_localized_vector(nonlin_sys,
138  *_base_sol).release());
139 
140  // also create localized solution vectos for the bassis vectors
141  std::vector<libMesh::NumericVector<Real>*> localized_basis(n_basis);
142  for (unsigned int i=0; i<n_basis; i++)
143  localized_basis[i] = build_localized_vector(nonlin_sys, *basis[i]).release();
144 
145 
146  // if a solution function is attached, initialize it
147  if (_sol_function && _base_sol)
149 
150 
151  libMesh::MeshBase::const_element_iterator el =
152  nonlin_sys.get_mesh().active_local_elements_begin();
153  const libMesh::MeshBase::const_element_iterator end_el =
154  nonlin_sys.get_mesh().active_local_elements_end();
155 
157  &ops = dynamic_cast<MAST::FluidStructureAssemblyElemOperations&>(*_elem_ops);
158 
159  for ( ; el != end_el; ++el) {
160 
161  const libMesh::Elem* elem = *el;
162 
163  dof_map.dof_indices (elem, dof_indices);
164 
165  // get the solution
166  unsigned int ndofs = (unsigned int)dof_indices.size();
167  sol.setZero(ndofs);
168  vec.setZero(ndofs);
169  mat.setZero(ndofs, ndofs);
170  basis_mat.setZero(ndofs, n_basis);
171 
172  for (unsigned int i=0; i<dof_indices.size(); i++) {
173 
174  if (_base_sol)
175  sol(i) = (*localized_solution)(dof_indices[i]);
176 
177 
178  for (unsigned int j=0; j<n_basis; j++)
179  basis_mat(i,j) = (*localized_basis[j])(dof_indices[i]);
180  }
181 
182 
183  // if (_sol_function)
184  // physics_elem->attach_active_solution_function(*_sol_function);
185 
186  MAST::GeomElem geom_elem;
187  _elem_ops->set_elem_data(elem->dim(), *elem, geom_elem);
188  geom_elem.init(*elem, *_system);
189 
190  _elem_ops->init(geom_elem);
192  _elem_ops->set_elem_velocity(vec); // set to zero value
193  _elem_ops->set_elem_acceleration(vec); // set to zero value
194 
195 
196  // now iterative over all qty types in the map and assemble them
197  it = mat_qty_map.begin();
198  end = mat_qty_map.end();
199 
200  for ( ; it != end; it++) {
201 
202  ops.set_qty_to_evaluate(it->first);
203  ops.elem_calculations(true, vec, mat);
204 
205  DenseRealMatrix m;
206  MAST::copy(m, mat);
207  dof_map.constrain_element_matrix(m, dof_indices);
208  MAST::copy(mat, m);
209 
210  // now add to the reduced order matrix
211  (*it->second) += basis_mat.transpose() * mat * basis_mat;
212  }
213 
215  // physics_elem->detach_active_solution_function();
216 
217  }
218 
219 
220  // if a solution function is attached, clear it
221  if (_sol_function)
222  _sol_function->clear();
223 
224 
225  // delete the localized basis vectors
226  for (unsigned int i=0; i<basis.size(); i++)
227  delete localized_basis[i];
228 
229  // sum the matrix and provide it to each processor
230  it = mat_qty_map.begin();
231  end = mat_qty_map.end();
232 
233 
234  for ( ; it != end; it++)
235  MAST::parallel_sum(_system->system().comm(), *(it->second));
236 }
237 
238 
239 
240 
241 
242 void
246  std::vector<libMesh::NumericVector<Real>*>& basis,
247  std::map<MAST::StructuralQuantityType, RealMatrixX*>& mat_qty_map) {
248 
249 
250  MAST::NonlinearSystem& nonlin_sys = _system->system();
251 
252  unsigned int
253  n_basis = (unsigned int)basis.size();
254 
255 
256  // initialize the quantities to zero matrices
257  std::map<MAST::StructuralQuantityType, RealMatrixX*>::iterator
258  it = mat_qty_map.begin(),
259  end = mat_qty_map.end();
260 
261  for ( ; it != end; it++)
262  *it->second = RealMatrixX::Zero(n_basis, n_basis);
263 
264  // iterate over each element, initialize it and get the relevant
265  // analysis quantities
266  RealVectorX vec, sol, dsol;
267  RealMatrixX mat, basis_mat;
268 
269  std::vector<libMesh::dof_id_type> dof_indices;
270  const libMesh::DofMap& dof_map = _system->system().get_dof_map();
271 
272 
273  std::unique_ptr<libMesh::NumericVector<Real> >
274  localized_solution,
275  localized_solution_sens;
276 
277  if (_base_sol) {
278 
279  // make sure that the solution sensitivity is provided
280  libmesh_assert(_base_sol_sensitivity);
281 
282  localized_solution.reset(build_localized_vector(nonlin_sys,
283  *_base_sol).release());
284  localized_solution_sens.reset(build_localized_vector(nonlin_sys,
285  *_base_sol_sensitivity).release());
286  }
287 
288  // also create localized solution vectos for the bassis vectors
289  std::vector<libMesh::NumericVector<Real>*> localized_basis(n_basis);
290  for (unsigned int i=0; i<n_basis; i++)
291  localized_basis[i] = build_localized_vector(nonlin_sys, *basis[i]).release();
292 
293 
294  // if a solution function is attached, initialize it
295  if (_sol_function && _base_sol)
297 
298 
299  libMesh::MeshBase::const_element_iterator el =
300  nonlin_sys.get_mesh().active_local_elements_begin();
301  const libMesh::MeshBase::const_element_iterator end_el =
302  nonlin_sys.get_mesh().active_local_elements_end();
303 
305  &ops = dynamic_cast<MAST::FluidStructureAssemblyElemOperations&>(*_elem_ops);
306 
307  for ( ; el != end_el; ++el) {
308 
309  const libMesh::Elem* elem = *el;
310 
311  dof_map.dof_indices (elem, dof_indices);
312 
313 
314  // get the solution
315  unsigned int ndofs = (unsigned int)dof_indices.size();
316  sol.setZero(ndofs);
317  dsol.setZero(ndofs);
318  vec.setZero(ndofs);
319  mat.setZero(ndofs, ndofs);
320  basis_mat.setZero(ndofs, n_basis);
321 
322  MAST::GeomElem geom_elem;
323  _elem_ops->set_elem_data(elem->dim(), *elem, geom_elem);
324  geom_elem.init(*elem, *_system);
325  _elem_ops->init(geom_elem);
326 
327  for (unsigned int i=0; i<dof_indices.size(); i++) {
328 
329  if (_base_sol) {
330 
332  sol(i) = (*localized_solution)(dof_indices[i]);
333  dsol(i) = (*localized_solution_sens)(dof_indices[i]);
334  }
335 
336  for (unsigned int j=0; j<n_basis; j++)
337  basis_mat(i,j) = (*localized_basis[j])(dof_indices[i]);
338  }
339 
342  _elem_ops->set_elem_velocity(vec); // set to zero value
343  _elem_ops->set_elem_acceleration(vec); // set to zero value
344 
345  // if (_sol_function)
346  // physics_elem->attach_active_solution_function(*_sol_function);
347 
348  // now iterative over all qty types in the map and assemble them
349  it = mat_qty_map.begin();
350  end = mat_qty_map.end();
351 
352  for ( ; it != end; it++) {
353 
354  ops.set_qty_to_evaluate(it->first);
355  ops.elem_sensitivity_calculations(f, true, vec, mat);
356 
357  DenseRealMatrix m;
358  MAST::copy(m, mat);
359  dof_map.constrain_element_matrix(m, dof_indices);
360  MAST::copy(mat, m);
361 
362  // now add to the reduced order matrix
363  (*it->second) += basis_mat.transpose() * mat * basis_mat;
364  }
365 
367  // physics_elem->detach_active_solution_function();
368 
369  }
370 
371 
372  // if a solution function is attached, clear it
373  if (_sol_function)
374  _sol_function->clear();
375 
376 
377  // delete the localized basis vectors
378  for (unsigned int i=0; i<basis.size(); i++)
379  delete localized_basis[i];
380 
381  // sum the matrix and provide it to each processor
382  it = mat_qty_map.begin();
383  end = mat_qty_map.end();
384 
385 
386  for ( ; it != end; it++)
387  MAST::parallel_sum(_system->system().comm(), *(it->second));
388 }
389 
MAST::AssemblyElemOperations * _elem_ops
provides assembly elem operations for use by this class
MAST::NonlinearSystem & system()
libMesh::DenseMatrix< Real > DenseRealMatrix
void clear_base_solution(bool if_sens=false)
Clears the pointer to base solution.
This class implements a system for solution of nonlinear systems.
virtual void set_elem_solution_sensitivity(const RealVectorX &sol)
sets the element solution sensitivity
virtual void assemble_reduced_order_quantity_sensitivity(const MAST::FunctionBase &f, std::vector< libMesh::NumericVector< Real > * > &basis, std::map< MAST::StructuralQuantityType, RealMatrixX * > &mat_qty_map)
calculates the sensitivity of reduced order matrix given the basis provided in basis.
virtual void set_elem_solution(const RealVectorX &sol)
sets the element solution
MAST::SystemInitialization * _system
System for which this assembly is performed.
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_discipline_and_system()
clears association with a system to this discipline
StructuralFluidInteractionAssembly()
constructor associates this assembly object with the system
Matrix< Real, Dynamic, Dynamic > RealMatrixX
void parallel_sum(const libMesh::Parallel::Communicator &c, RealMatrixX &mat)
Definition: utility.h:214
virtual void assemble_reduced_order_quantity(std::vector< libMesh::NumericVector< Real > * > &basis, std::map< MAST::StructuralQuantityType, RealMatrixX * > &mat_qty_map)
calculates the reduced order matrix given the basis provided in basis.
virtual ~StructuralFluidInteractionAssembly()
destructor resets the association of this assembly object with the system
virtual void clear_discipline_and_system()
clears association with a system to this discipline, and vice-a-versa
void copy(DenseRealMatrix &m1, const RealMatrixX &m2)
Definition: utility.h:167
void clear()
clear the solution
Matrix< Real, Dynamic, 1 > RealVectorX
void set_base_solution(const libMesh::NumericVector< Real > &sol, bool if_sens=false)
if the eigenproblem is defined about a non-zero base solution, then this method provides the object w...
virtual void init(const MAST::GeomElem &elem)=0
initializes the object for calculation of element quantities for the specified elem.
virtual void elem_calculations(bool if_jac, RealVectorX &vec, RealMatrixX &mat)
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...
This class acts as a wrapper around libMesh::Elem for the purpose of providing a uniform interface fo...
Definition: geom_elem.h:59
virtual void set_elem_velocity(const RealVectorX &vel)
sets the element velocity
const libMesh::NumericVector< Real > * _base_sol_sensitivity
sensitivity of base solution may be needed for sensitivity analysis.
virtual void elem_sensitivity_calculations(const MAST::FunctionBase &f, bool if_jac, RealVectorX &vec, RealMatrixX &mat)
void init(const libMesh::NumericVector< Real > &sol, const libMesh::NumericVector< Real > *dsol=nullptr)
initializes the data structures to perform the interpolation function of sol.
const libMesh::NumericVector< Real > * _base_sol
base solution about which this eigenproblem is defined.
virtual void clear_elem()
clears the element initialization
void use_base_sol_for_sensitivity(bool f)
if set to true, the sensitivity calculation will include the sensitivity of base solution.
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 void set_elem_acceleration(const RealVectorX &accel)
sets the element acceleration
MAST::MeshFieldFunction * _sol_function
system solution that will be initialized before each solution