libMesh
quadrature_composite.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 #include "libmesh/libmesh_config.h"
20 #if defined(LIBMESH_HAVE_TRIANGLE) && defined(LIBMESH_HAVE_TETGEN)
21 
22 #include "libmesh/fe.h"
23 #include "libmesh/quadrature_gauss.h"
24 #include "libmesh/quadrature_trap.h"
25 #include "libmesh/quadrature_simpson.h"
26 #include "libmesh/quadrature_composite.h"
27 #include "libmesh/elem.h"
28 
29 
30 
31 namespace libMesh
32 {
33 
34 
35 template <class QSubCell>
36 QComposite<QSubCell>::QComposite(const unsigned int d,
37  const Order o) :
38  QSubCell(d,o), // explicitly call base class constructor
39  _q_subcell(d,o),
40  _lagrange_fe(FEBase::build (d, FEType (FIRST, LAGRANGE)))
41 {
42  // explicitly call the init function in 1D since the
43  // other tensor-product rules require this one.
44  // note that EDGE will not be used internally, however
45  // if we called the function with INVALID_ELEM it would try to
46  // be smart and return, thinking it had already done the work.
47  if (_dim == 1)
49 
50  libmesh_assert (_lagrange_fe.get() != libmesh_nullptr);
51 
52  _lagrange_fe->attach_quadrature_rule (&_q_subcell);
53 }
54 
55 
56 
57 template <class QSubCell>
58 QComposite<QSubCell>::~QComposite()
59 {}
60 
61 
62 
63 template <class QSubCell>
64 void QComposite<QSubCell>::init (const Elem & elem,
65  const std::vector<Real> & vertex_distance_func,
66  unsigned int p_level)
67 {
68  libmesh_assert_equal_to (vertex_distance_func.size(), elem.n_vertices());
69  libmesh_assert_equal_to (_dim, elem.dim());
70 
71  // if we are not cut, revert to simple base class init() method.
72  if (!_elem_cutter.is_cut (elem, vertex_distance_func))
73  {
74  _q_subcell.init (elem.type(), p_level);
75  _points = _q_subcell.get_points();
76  _weights = _q_subcell.get_weights();
77 
78  //this->print_info();
79  return;
80  }
81 
82  // Get a pointer to the element's reference element. We want to
83  // perform cutting on the reference element such that the quadrature
84  // point locations of the subelements live in the reference
85  // coordinate system, thereby eliminating the need for inverse
86  // mapping.
87  const Elem * reference_elem = elem.reference_elem();
88 
89  libmesh_assert (reference_elem != libmesh_nullptr);
90 
91  _elem_cutter(*reference_elem, vertex_distance_func);
92  //_elem_cutter(elem, vertex_distance_func);
93 
94  // clear our state & accumulate points from subelements
95  _points.clear();
96  _weights.clear();
97 
98  // inside subelem
99  {
100  const std::vector<Elem const *> & inside_elem (_elem_cutter.inside_elements());
101  std::cout << inside_elem.size() << " elements inside\n";
102 
103  this->add_subelem_values(inside_elem);
104  }
105 
106  // outside subelem
107  {
108  const std::vector<Elem const *> & outside_elem (_elem_cutter.outside_elements());
109  std::cout << outside_elem.size() << " elements outside\n";
110 
111  this->add_subelem_values(outside_elem);
112  }
113 
114  this->print_info();
115 }
116 
117 
118 
119 template <class QSubCell>
120 void QComposite<QSubCell>::add_subelem_values (const std::vector<Elem const *> & subelem)
121 
122 {
123  const std::vector<Real> & subelem_weights = _lagrange_fe->get_JxW();
124  const std::vector<Point> & subelem_points = _lagrange_fe->get_xyz();
125 
126  for (std::vector<Elem const *>::const_iterator it = subelem.begin();
127  it!=subelem.end(); ++it)
128  {
129  // tetgen seems to create 0-volume cells on occasion, but we *should*
130  // be catching that appropriately now inside the ElemCutter class.
131  // Just in case trap here, describe the error, and abort.
132  libmesh_try
133  {
134  _lagrange_fe->reinit(*it);
135  _weights.insert(_weights.end(),
136  subelem_weights.begin(), subelem_weights.end());
137 
138  _points.insert(_points.end(),
139  subelem_points.begin(), subelem_points.end());
140  }
141  libmesh_catch (...)
142  {
143  libMesh::err << "ERROR: found a bad cut cell!\n";
144 
145  for (unsigned int n=0; n<(*it)->n_nodes(); n++)
146  libMesh::err << (*it)->point(n) << std::endl;
147 
148  libmesh_error_msg("Tetgen may have created a 0-volume cell during Cutcell integration.");
149  }
150  }
151 }
152 
153 
154 
155 //--------------------------------------------------------------
156 // Explicit instantiations
157 template class QComposite<QGauss>;
158 template class QComposite<QTrap>;
159 template class QComposite<QSimpson>;
160 
161 } // namespace libMesh
162 
163 #endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_HAVE_TETGEN
OStreamProxy err
const unsigned int _dim
The spatial dimension of the quadrature rule.
Definition: quadrature.h:311
const class libmesh_nullptr_t libmesh_nullptr
The libMesh namespace provides an interface to certain functionality in the library.
std::vector< Point > _points
The locations of the quadrature points in reference element space.
Definition: quadrature.h:335
std::vector< Real > _weights
The quadrature weights.
Definition: quadrature.h:341
libmesh_assert(j)
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
void print_info(std::ostream &os=libMesh::out) const
Prints information relevant to the quadrature rule, by default to libMesh::out.
Definition: quadrature.h:364
FEGenericBase< Real > FEBase
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32