Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
20 : // C++ includes
21 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 : #include <cmath> // for std::sqrt
23 :
24 :
25 : // Local includes
26 : #include "libmesh/libmesh_common.h"
27 : #include "libmesh/fe.h"
28 : #include "libmesh/quadrature.h"
29 : #include "libmesh/elem.h"
30 : #include "libmesh/libmesh_logging.h"
31 :
32 : namespace libMesh
33 : {
34 :
35 :
36 :
37 :
38 : //-------------------------------------------------------
39 : // Full specialization for 0D when this is a useless method
40 : template <>
41 0 : void FEXYZ<0>::reinit(const Elem *,
42 : const unsigned int,
43 : const Real,
44 : const std::vector<Point> * const,
45 : const std::vector<Real> * const)
46 : {
47 0 : libmesh_error_msg("ERROR: This method only makes sense for 2D/3D elements!");
48 : }
49 :
50 :
51 : //-------------------------------------------------------
52 : // Full specialization for 1D when this is a useless method
53 : template <>
54 0 : void FEXYZ<1>::reinit(const Elem *,
55 : const unsigned int,
56 : const Real,
57 : const std::vector<Point> * const,
58 : const std::vector<Real> * const)
59 : {
60 0 : libmesh_error_msg("ERROR: This method only makes sense for 2D/3D elements!");
61 : }
62 :
63 :
64 : //-------------------------------------------------------
65 : // Method for 2D, 3D
66 : template <unsigned int Dim>
67 1721016 : void FEXYZ<Dim>::reinit(const Elem * elem,
68 : const unsigned int s,
69 : const Real,
70 : const std::vector<Point> * const pts,
71 : const std::vector<Real> * const weights)
72 : {
73 143418 : libmesh_assert(elem);
74 143418 : libmesh_assert (this->qrule != nullptr || pts != nullptr);
75 : // We don't do this for 1D elements!
76 : libmesh_assert_not_equal_to (Dim, 1);
77 :
78 : // Build the side of interest
79 1721016 : const std::unique_ptr<const Elem> side(elem->build_side_ptr(s));
80 :
81 : // Find the max p_level to select
82 : // the right quadrature rule for side integration
83 1721016 : unsigned int side_p_level = elem->p_level();
84 1864434 : if (elem->neighbor_ptr(s) != nullptr)
85 1776086 : side_p_level = std::max(side_p_level, elem->neighbor_ptr(s)->p_level());
86 :
87 : // Initialize the shape functions at the user-specified
88 : // points
89 1721016 : if (pts != nullptr)
90 : {
91 : // We can't get away without recomputing shape functions next
92 : // time
93 0 : this->shapes_on_quadrature = false;
94 :
95 : // Set the element type
96 0 : this->_elem = elem;
97 0 : this->_elem_type = elem->type();
98 :
99 : // Initialize the face shape functions
100 0 : this->_fe_map->template init_face_shape_functions<Dim>(*pts, side.get());
101 0 : if (weights != nullptr)
102 : {
103 0 : this->compute_face_values (elem, side.get(), *weights);
104 : }
105 : else
106 : {
107 0 : std::vector<Real> dummy_weights (pts->size(), 1.);
108 : // Compute data on the face for integration
109 0 : this->compute_face_values (elem, side.get(), dummy_weights);
110 : }
111 : }
112 : else
113 : {
114 : // initialize quadrature rule
115 1721016 : this->qrule->init(*side, side_p_level);
116 :
117 : {
118 : // Set the element
119 1721016 : this->_elem = elem;
120 1721016 : this->_elem_type = elem->type();
121 :
122 : // Initialize the face shape functions
123 1864434 : this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(), side.get());
124 : }
125 : // We can't get away without recomputing shape functions next
126 : // time
127 1721016 : this->shapes_on_quadrature = false;
128 : // Compute data on the face for integration
129 1864434 : this->compute_face_values (elem, side.get(), this->qrule->get_weights());
130 : }
131 1721016 : }
132 :
133 :
134 :
135 : template <unsigned int Dim>
136 1721016 : void FEXYZ<Dim>::compute_face_values(const Elem * elem,
137 : const Elem * side,
138 : const std::vector<Real> & qw)
139 : {
140 143418 : libmesh_assert(elem);
141 143418 : libmesh_assert(side);
142 :
143 286836 : LOG_SCOPE("compute_face_values()", "FEXYZ");
144 :
145 : // The number of quadrature points.
146 286836 : const std::size_t n_qp = qw.size();
147 :
148 : // Number of shape functions in the finite element approximation
149 : // space.
150 286836 : const unsigned int n_approx_shape_functions =
151 1434180 : this->n_dofs(elem, this->get_order());
152 :
153 : // Resize the shape functions and their gradients
154 1721016 : this->phi.resize (n_approx_shape_functions);
155 1721016 : this->dphi.resize (n_approx_shape_functions);
156 1721016 : this->dphidx.resize (n_approx_shape_functions);
157 1721016 : this->dphidy.resize (n_approx_shape_functions);
158 1721016 : this->dphidz.resize (n_approx_shape_functions);
159 :
160 9349776 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
161 : {
162 8264490 : this->phi[i].resize (n_qp);
163 8264490 : this->dphi[i].resize (n_qp);
164 8264490 : this->dphidx[i].resize (n_qp);
165 8264490 : this->dphidy[i].resize (n_qp);
166 8264490 : this->dphidz[i].resize (n_qp);
167 : }
168 :
169 1721016 : this->_fe_map->compute_face_map(this->dim, qw, side);
170 :
171 143418 : const std::vector<libMesh::Point> & xyz = this->_fe_map->get_xyz();
172 :
173 1721016 : switch (this->dim)
174 : {
175 : // A 2D finite element living in either 2D or 3D space.
176 : // This means the boundary is a 1D finite element, i.e.
177 : // and EDGE2 or EDGE3.
178 19770 : case 2:
179 : {
180 : // compute the shape function values & gradients
181 1046160 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
182 2426760 : for (std::size_t p=0; p<n_qp; p++)
183 : {
184 1752660 : this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
185 :
186 1752660 : this->dphi[i][p](0) =
187 1752660 : this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
188 :
189 1752660 : this->dphi[i][p](1) =
190 1752660 : this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
191 :
192 : #if LIBMESH_DIM == 3
193 1617840 : this->dphi[i][p](2) = // can only assign to the Z component if LIBMESH_DIM==3
194 : #endif
195 1752660 : this->dphidz[i][p] = 0.;
196 : }
197 :
198 : // done computing face values
199 19770 : break;
200 : }
201 :
202 : // A 3D finite element living in 3D space.
203 123648 : case 3:
204 : {
205 : // compute the shape function values & gradients
206 8303616 : for (unsigned int i=0; i<n_approx_shape_functions; i++)
207 34099200 : for (std::size_t p=0; p<n_qp; p++)
208 : {
209 29552640 : this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
210 :
211 29552640 : this->dphi[i][p](0) =
212 29552640 : this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
213 :
214 29552640 : this->dphi[i][p](1) =
215 29552640 : this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
216 :
217 31825920 : this->dphi[i][p](2) =
218 29552640 : this->dphidz[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 2, xyz[p]);
219 : }
220 :
221 : // done computing face values
222 123648 : break;
223 : }
224 :
225 0 : default:
226 0 : libmesh_error_msg("Invalid dim " << this->dim);
227 : }
228 1721016 : }
229 :
230 :
231 :
232 :
233 : //--------------------------------------------------------------
234 : // Explicit instantiations (doesn't make sense in 1D!) using fe_macro.h's macro
235 :
236 : template class LIBMESH_EXPORT FEXYZ<2>;
237 : template class LIBMESH_EXPORT FEXYZ<3>;
238 :
239 :
240 : } // namespace libMesh
|
