GRINS-0.8.0
parsed_boundary_qoi.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2017 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 
26 // This class
28 
29 // GRINS
30 #include "grins/multiphysics_sys.h"
31 #include "grins/assembly_context.h"
32 
33 // libMesh
34 #include "libmesh/getpot.h"
35 #include "libmesh/fem_system.h"
36 #include "libmesh/quadrature.h"
37 #include "libmesh/fe_base.h"
38 #include "libmesh/parsed_fem_function.h"
39 
40 namespace GRINS
41 {
42  ParsedBoundaryQoI::ParsedBoundaryQoI( const std::string& qoi_name )
43  : QoIBase(qoi_name) {}
44 
46  : QoIBase(original)
47  {
48  if (original.qoi_functional.get())
49  {
50  this->qoi_functional = original.qoi_functional->clone();
51  this->move_parameter
53  (original.qoi_functional.get()),
55  (this->qoi_functional.get()));
56  }
57 
58  this->_bc_ids = original._bc_ids;
59  }
60 
62 
64  {
65  return new ParsedBoundaryQoI( *this );
66  }
67 
68  void ParsedBoundaryQoI::init( const GetPot& input,
69  const MultiphysicsSystem& system,
70  unsigned int /*qoi_num*/ )
71  {
72  // Read boundary ids on which we want to compute qoi
73  int num_bcs = input.vector_variable_size("QoI/ParsedBoundary/bc_ids");
74 
75  if( num_bcs <= 0 )
76  {
77  std::cerr << "Error: Must specify at least one boundary id to compute"
78  << " parsed boundary QoI." << std::endl
79  << "Found: " << num_bcs << std::endl;
80  libmesh_error();
81  }
82 
83  for( int i = 0; i < num_bcs; i++ )
84  {
85  _bc_ids.insert( input("QoI/ParsedBoundary/bc_ids", -1, i ) );
86  }
87 
90  (system, ""));
91  this->qoi_functional.reset(qf);
92 
93  this->set_parameter(*qf, input,
94  "QoI/ParsedBoundary/qoi_functional", "DIE!");
95  }
96 
98  {
99  libMesh::FEBase* side_fe;
100  context.get_side_fe<libMesh::Real>(0, side_fe);
101  side_fe->get_JxW();
102  side_fe->get_xyz();
103 
104  qoi_functional->init_context(context);
105  }
106 
108  const unsigned int qoi_index )
109  {
110  bool on_correct_side = false;
111 
112  for (std::set<libMesh::boundary_id_type>::const_iterator id =
113  _bc_ids.begin(); id != _bc_ids.end(); id++ )
114  if( context.has_side_boundary_id( (*id) ) )
115  {
116  on_correct_side = true;
117  break;
118  }
119 
120  if (!on_correct_side)
121  return;
122 
123  libMesh::FEBase* side_fe;
124  context.get_side_fe<libMesh::Real>(0, side_fe);
125  const std::vector<libMesh::Real> &JxW = side_fe->get_JxW();
126 
127  const std::vector<libMesh::Point>& x_qp = side_fe->get_xyz();
128 
129  unsigned int n_qpoints = context.get_side_qrule().n_points();
130 
132  libMesh::Number& qoi = context.get_qois()[qoi_index];
133 
134  for( unsigned int qp = 0; qp != n_qpoints; qp++ )
135  {
136  const libMesh::Number func_val =
137  (*qoi_functional)(context, x_qp[qp], context.get_time());
138 
139  qoi += func_val * JxW[qp];
140  }
141  }
142 
144  const unsigned int qoi_index )
145  {
146  bool on_correct_side = false;
147 
148  for (std::set<libMesh::boundary_id_type>::const_iterator id =
149  _bc_ids.begin(); id != _bc_ids.end(); id++ )
150  if( context.has_side_boundary_id( (*id) ) )
151  {
152  on_correct_side = true;
153  break;
154  }
155 
156  if (!on_correct_side)
157  return;
158 
159  libMesh::FEBase* side_fe;
160  context.get_side_fe<libMesh::Real>(0, side_fe);
161  const std::vector<libMesh::Real> &JxW = side_fe->get_JxW();
162 
163  const std::vector<libMesh::Point>& x_qp = side_fe->get_xyz();
164 
165  // Local DOF count and quadrature point count
166  const unsigned int n_u_dofs = context.get_dof_indices().size();
167 
168  unsigned int n_qpoints = context.get_side_qrule().n_points();
169 
170  // Local solution vector - non-const version for finite
171  // differenting purposes
172  libMesh::DenseVector<libMesh::Number>& elem_solution =
173  const_cast<libMesh::DenseVector<libMesh::Number>&>
174  (context.get_elem_solution());
175 
177  libMesh::DenseVector<libMesh::Number> &Qu =
178  context.get_qoi_derivatives()[qoi_index];
179 
180  for( unsigned int qp = 0; qp != n_qpoints; qp++ )
181  {
182  // Central finite differencing to approximate derivatives.
183  // FIXME - we should hook the FParserAD stuff into
184  // ParsedFEMFunction
185 
186  for( unsigned int i = 0; i != n_u_dofs; ++i )
187  {
188  libMesh::Number &current_solution = elem_solution(i);
189  const libMesh::Number original_solution = current_solution;
190 
191  current_solution = original_solution + libMesh::TOLERANCE;
192 
193  const libMesh::Number plus_val =
194  (*qoi_functional)(context, x_qp[qp], context.get_time());
195 
196  current_solution = original_solution - libMesh::TOLERANCE;
197 
198  const libMesh::Number minus_val =
199  (*qoi_functional)(context, x_qp[qp], context.get_time());
200 
201  Qu(i) += (plus_val - minus_val) *
202  (0.5 / libMesh::TOLERANCE) * JxW[qp];
203 
204  // Don't forget to restore the correct solution...
205  current_solution = original_solution;
206  }
207  }
208  }
209 
210 } //namespace GRINS
virtual void set_parameter(libMesh::Number &param_variable, const GetPot &input, const std::string &param_name, libMesh::Number param_default)
Each subclass can simultaneously read a parameter value from.
virtual void init_context(AssemblyContext &context)
ParsedBoundaryQoI()
User never call default constructor.
virtual QoIBase * clone() const
Required to provide clone (deep-copy) for adding QoI object to libMesh objects.
GRINS namespace.
Parsed Boundary QoI.
std::set< libMesh::boundary_id_type > _bc_ids
List of boundary ids on which we want to compute this QoI.
virtual void init(const GetPot &input, const MultiphysicsSystem &system, unsigned int qoi_num)
Initialize local variables.
virtual void side_qoi_derivative(AssemblyContext &context, const unsigned int qoi_index)
Compute the qoi derivative with respect to the solution.
Interface with libMesh for solving Multiphysics problems.
libMesh::UniquePtr< libMesh::FEMFunctionBase< libMesh::Number > > qoi_functional
virtual void side_qoi(AssemblyContext &context, const unsigned int qoi_index)
Compute the qoi value.
virtual void move_parameter(const libMesh::Number &old_parameter, libMesh::Number &new_parameter)
When cloning an object, we need to update parameter pointers.

Generated on Tue Dec 19 2017 12:47:28 for GRINS-0.8.0 by  doxygen 1.8.9.1