GRINS-0.8.0
distance_function.h
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 #ifndef GRINS_DISTANCE_FUNCTION_H
26 #define GRINS_DISTANCE_FUNCTION_H
27 
28 // system
29 //#include <limits>
30 
31 // libmesh
32 #include "libmesh/libmesh_base.h"
33 #include "libmesh/libmesh_common.h"
34 #include "libmesh/auto_ptr.h"
35 #include "libmesh/dense_vector.h"
36 #include "libmesh/dense_matrix.h"
37 #include "libmesh/fe_base.h"
38 #include "libmesh/system.h"
39 #include "libmesh/elem.h"
40 
41 // Forward Declarations
42 namespace libMesh {
43  class EquationSystems;
44  class BoundaryMesh;
45  class Node;
46 }
47 
48 
49 namespace GRINS {
50 
51  // This class provides the functionality to compute the distance to
52  // the nearest no slip wall boundary.
53  class DistanceFunction : public libMesh::System::Initialization
54  {
55  public:
56 
60  DistanceFunction (libMesh::EquationSystems &es_in, const libMesh::UnstructuredMesh &bm_in);
61 
66 
70  virtual void initialize ();
71 
75  libMesh::Real node_to_boundary (const libMesh::Node* node);
76 
81  void compute ();
82 
86  libMesh::UniquePtr< libMesh::DenseVector<libMesh::Real> > interpolate (const libMesh::Elem* elem, const std::vector<libMesh::Point>& qts) const;
87 
88 
89 
90  private:
91 
95  libMesh::EquationSystems &_equation_systems;
96 
100  const libMesh::UnstructuredMesh &_boundary_mesh;
101 
110  libMesh::UniquePtr<libMesh::FEBase> _dist_fe;
111 
112  };
113 
119  {
120  public:
121 
122  // ctor
124 
125  // dtor
127 
131  template <class f>
132  void operator()(const libMesh::DenseVector<libMesh::Real> &U,
133  const libMesh::DenseVector<libMesh::Real> & F0,
134  f &F,
135  libMesh::DenseMatrix<libMesh::Real> &dFdU)
136  {
137 
138  Up = U;
139  Fp.resize(U.size());
140 
141  dFdU.resize(U.size(),U.size());
142 
143  // compute dF_i/dU_j.
144  for (unsigned int j=0; j<U.size(); j++)
145  {
146  // tell the user's function which component we are perturbing
147  // at this iteration. this information may be used to avoid
148  // repeated calculations of values which are known to be
149  // unaffected for this perturbation.
150  //F.set_perturbed_component(j);
151 
152  // define the pertubation for this component
153  const libMesh::Real
154  pert = 4.e-8*(std::max (std::abs(U(j)), 1.e-3)); // U can be 0, pertubation cannot
155 
156  Up(j) += pert;
157 
158  const libMesh::Real // note this difference may not strictly be pert
159  invpert = 1./(Up(j) - U(j)); // due to truncation error in the preceeding sum.
160 
161  // evaluate F at the perturbed state
162  F(Up,Fp);
163 
164  // remove perturbation for next evaluation
165  Up(j) = U(j);
166 
167  for (unsigned int i=0; i<U.size(); i++)
168  dFdU(i,j) = (Fp(i) - F0(i))*invpert;
169  }
170  }
171 
172  private:
173 
174  // work vectors
175  libMesh::DenseVector<libMesh::Real> Up, Fp, Um, Fm;
176  };
177 
178 
179 
189  {
190  public:
191 
192  // ctor
193  ComputeDistanceResidual(const libMesh::Elem* belem, const libMesh::Point* point) :
194  _belem(*belem),
195  _dim(_belem.dim()+1),
196  _p(*point),
197  _fe (libMesh::FEBase::build(_belem.dim(), libMesh::FEType(libMesh::FIRST, libMesh::LAGRANGE))),
198  fe (_fe.get())
199  {
200  // only supporting QUAD4 this way for now
201  libmesh_assert(_belem.type()==libMesh::QUAD4);
202  }
203 
204  // dtor
206 
207  // Calculate the residual
208  void operator()(const libMesh::DenseVector<libMesh::Real> &U,
209  libMesh::DenseVector<libMesh::Real> &F)
210  {
211 
212  libmesh_assert(U.size()==_dim-1);
213  libmesh_assert(F.size()==_dim-1);
214 
215  libMesh::DenseVector<libMesh::Real> xx(_dim);
216  xx.zero();
217 
218  libMesh::DenseMatrix<libMesh::Real> xx_U(_dim, _dim-1);
219  xx_U.zero();
220 
221  std::vector<libMesh::Point> xi(1);
222 
223  xi[0](0) = U(0);
224  xi[0](1) = U(1);
225 
226  // grab basis functions (evaluated at qpts)
227  const std::vector<std::vector<libMesh::Real> > &basis = fe->get_phi();
228  const std::vector<std::vector<libMesh::RealGradient> > &dbasis = fe->get_dphi();
229 
230  // reinitialize finite element data at xi
231  fe->reinit(&_belem, &xi);
232 
233  // interpolate location and derivatives
234  for (unsigned int inode=0; inode<_belem.n_nodes(); ++inode)
235  for (unsigned int idim=0; idim<_dim; ++idim)
236  {
237  xx(idim) += _belem.point(inode)(idim) * basis[inode][0];
238 
239  for (unsigned int jdim=0; jdim<_dim-1; ++jdim)
240  xx_U(idim, jdim) += _belem.point(inode)(idim) * dbasis[inode][0](jdim);
241  }
242 
243 
244  // form the residual
245  F.zero();
246  for (unsigned int ires=0; ires<_dim-1; ++ires)
247  for (unsigned int idim=0; idim<_dim; ++idim)
248  F(ires) += ( xx(idim) - _p(idim) ) * xx_U(idim, ires);
249 
250 
251  }
252 
253  private:
254 
255  // Reference to boundary element
256  const libMesh::Elem& _belem;
257  const unsigned int _dim;
258  const libMesh::Point& _p;
259 
260  libMesh::UniquePtr<libMesh::FEBase> _fe;
261  libMesh::FEBase *fe;
262 
263  // work vectors
264  libMesh::DenseVector<libMesh::Real> Up, Fp, Um, Fm;
265  };
266 
267 } // end namespace GRINS
268 
269 #endif // __distance_function_h__
libMesh::UniquePtr< libMesh::FEBase > _dist_fe
libMesh::DenseVector< libMesh::Real > Up
libMesh::DenseVector< libMesh::Real > Um
libMesh::EquationSystems & _equation_systems
DistanceFunction(libMesh::EquationSystems &es_in, const libMesh::UnstructuredMesh &bm_in)
GRINS namespace.
libMesh::DenseVector< libMesh::Real > Fp
void operator()(const libMesh::DenseVector< libMesh::Real > &U, const libMesh::DenseVector< libMesh::Real > &F0, f &F, libMesh::DenseMatrix< libMesh::Real > &dFdU)
ComputeDistanceResidual(const libMesh::Elem *belem, const libMesh::Point *point)
libMesh::UniquePtr< libMesh::DenseVector< libMesh::Real > > interpolate(const libMesh::Elem *elem, const std::vector< libMesh::Point > &qts) const
libMesh::UniquePtr< libMesh::FEBase > _fe
libMesh::DenseVector< libMesh::Real > Fm
libMesh::DenseVector< libMesh::Real > Fp
libMesh::DenseVector< libMesh::Real > Um
libMesh::DenseVector< libMesh::Real > Fm
libMesh::Real node_to_boundary(const libMesh::Node *node)
void operator()(const libMesh::DenseVector< libMesh::Real > &U, libMesh::DenseVector< libMesh::Real > &F)
const libMesh::UnstructuredMesh & _boundary_mesh
libMesh::DenseVector< libMesh::Real > Up

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