Intrepid2
Intrepid2_HGRAD_TET_Cn_FEM.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
49#ifndef __INTREPID2_HGRAD_TET_CN_FEM_HPP__
50#define __INTREPID2_HGRAD_TET_CN_FEM_HPP__
51
52#include "Intrepid2_Basis.hpp"
55
57#include "Teuchos_LAPACK.hpp"
58
59
60namespace Intrepid2 {
61
87 namespace Impl {
88
93 public:
94 typedef struct Tetrahedron<4> cell_topology_type;
98 template<EOperator opType>
99 struct Serial {
100 template<typename outputValueViewType,
101 typename inputPointViewType,
102 typename workViewType,
103 typename vinvViewType>
105 static void
109 const vinvViewType vinv );
110 };
111
112 template<typename DeviceType, ordinal_type numPtsPerEval,
113 typename outputValueValueType, class ...outputValueProperties,
114 typename inputPointValueType, class ...inputPointProperties,
115 typename vinvValueType, class ...vinvProperties>
116 static void
117 getValues( Kokkos::DynRankView<outputValueValueType,outputValueProperties...> outputValues,
118 const Kokkos::DynRankView<inputPointValueType, inputPointProperties...> inputPoints,
119 const Kokkos::DynRankView<vinvValueType, vinvProperties...> vinv,
120 const EOperator operatorType);
121
125 template<typename outputValueViewType,
126 typename inputPointViewType,
127 typename vinvViewType,
128 typename workViewType,
129 EOperator opType,
130 ordinal_type numPtsEval>
131 struct Functor {
132 outputValueViewType _outputValues;
133 const inputPointViewType _inputPoints;
134 const vinvViewType _vinv;
135 workViewType _work;
136
139 inputPointViewType inputPoints_,
140 vinvViewType vinv_,
142 : _outputValues(outputValues_), _inputPoints(inputPoints_),
143 _vinv(vinv_), _work(work_) {}
144
146 void operator()(const size_type iter) const {
147 const auto ptBegin = Util<ordinal_type>::min(iter*numPtsEval, _inputPoints.extent(0));
148 const auto ptEnd = Util<ordinal_type>::min(ptBegin+numPtsEval, _inputPoints.extent(0));
149
150 const auto ptRange = Kokkos::pair<ordinal_type,ordinal_type>(ptBegin, ptEnd);
151 const auto input = Kokkos::subview( _inputPoints, ptRange, Kokkos::ALL() );
152
153 typename workViewType::pointer_type ptr = _work.data() + _work.extent(0)*ptBegin*get_dimension_scalar(_work);
154
155 auto vcprop = Kokkos::common_view_alloc_prop(_work);
156 workViewType work(Kokkos::view_wrap(ptr,vcprop), (ptEnd-ptBegin)*_work.extent(0));
157
158 switch (opType) {
159 case OPERATOR_VALUE : {
160 auto output = Kokkos::subview( _outputValues, Kokkos::ALL(), ptRange );
161 Serial<opType>::getValues( output, input, work, _vinv );
162 break;
163 }
164 case OPERATOR_GRAD :
165 case OPERATOR_D1 :
166 case OPERATOR_D2 :
167 //case OPERATOR_D3 :
168 {
169 auto output = Kokkos::subview( _outputValues, Kokkos::ALL(), ptRange, Kokkos::ALL() );
170 Serial<opType>::getValues( output, input, work, _vinv );
171 break;
172 }
173 default: {
174 INTREPID2_TEST_FOR_ABORT( true,
175 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM::Functor) operator is not supported");
176
177 }
178 }
179 }
180 };
181 };
182 }
183
184 template<typename DeviceType = void,
185 typename outputValueType = double,
186 typename pointValueType = double>
188 : public Basis<DeviceType,outputValueType,pointValueType> {
189 public:
193
197
199
200 private:
201
204 Kokkos::DynRankView<scalarType,DeviceType> vinv_;
205
207 EPointType pointType_;
208
209 public:
210
213 Basis_HGRAD_TET_Cn_FEM(const ordinal_type order,
214 const EPointType pointType = POINTTYPE_EQUISPACED);
215
216
217
218 using Basis<DeviceType,outputValueType,pointValueType>::getValues;
219
220 virtual
221 void
222 getValues( OutputViewType outputValues,
223 const PointViewType inputPoints,
224 const EOperator operatorType = OPERATOR_VALUE) const override {
225#ifdef HAVE_INTREPID2_DEBUG
227 inputPoints,
228 operatorType,
229 this->getBaseCellTopology(),
230 this->getCardinality() );
231#endif
232 constexpr ordinal_type numPtsPerEval = Parameters::MaxNumPtsPerBasisEval;
233 Impl::Basis_HGRAD_TET_Cn_FEM::
234 getValues<DeviceType,numPtsPerEval>( outputValues,
235 inputPoints,
236 this->vinv_,
237 operatorType);
238 }
239
240 virtual
241 void
242 getDofCoords( ScalarViewType dofCoords ) const override {
243#ifdef HAVE_INTREPID2_DEBUG
244 // Verify rank of output array.
245 INTREPID2_TEST_FOR_EXCEPTION( dofCoords.rank() != 2, std::invalid_argument,
246 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM::getDofCoords) rank = 2 required for dofCoords array");
247 // Verify 0th dimension of output array.
248 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(dofCoords.extent(0)) != this->getCardinality(), std::invalid_argument,
249 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM::getDofCoords) mismatch in number of dof and 0th dimension of dofCoords array");
250 // Verify 1st dimension of output array.
251 INTREPID2_TEST_FOR_EXCEPTION( dofCoords.extent(1) != this->getBaseCellTopology().getDimension(), std::invalid_argument,
252 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM::getDofCoords) incorrect reference cell (1st) dimension in dofCoords array");
253#endif
254 Kokkos::deep_copy(dofCoords, this->dofCoords_);
255 }
256
257 virtual
258 void
259 getDofCoeffs( ScalarViewType dofCoeffs ) const override {
260#ifdef HAVE_INTREPID2_DEBUG
261 // Verify rank of output array.
262 INTREPID2_TEST_FOR_EXCEPTION( dofCoeffs.rank() != 1, std::invalid_argument,
263 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM::getdofCoeffs) rank = 1 required for dofCoeffs array");
264 // Verify 0th dimension of output array.
265 INTREPID2_TEST_FOR_EXCEPTION( static_cast<ordinal_type>(dofCoeffs.extent(0)) != this->getCardinality(), std::invalid_argument,
266 ">>> ERROR: (Intrepid2::Basis_HGRAD_TET_Cn_FEM::getdofCoeffs) mismatch in number of dof and 0th dimension of dofCoeffs array");
267#endif
268 Kokkos::deep_copy(dofCoeffs, 1.0);
269 }
270
271
272 void
273 getVandermondeInverse( ScalarViewType vinv ) const {
274 // has to be same rank and dimensions
275 Kokkos::deep_copy(vinv, this->vinv_);
276 }
277
278 virtual
279 const char*
280 getName() const override {
281 return "Intrepid2_HGRAD_TET_Cn_FEM";
282 }
283
284 virtual
285 bool
286 requireOrientation() const override {
287 return (this->basisDegree_ > 2);
288 }
289
290 Kokkos::DynRankView<typename ScalarViewType::const_value_type,DeviceType>
291
292 getVandermondeInverse() const {
293 return vinv_;
294 }
295
296 ordinal_type
297 getWorkSizePerPoint(const EOperator operatorType) const {
298 auto cardinality = getPnCardinality<3>(this->basisDegree_);
299 switch (operatorType) {
300 case OPERATOR_GRAD:
301 case OPERATOR_CURL:
302 case OPERATOR_D1:
303 return 7*cardinality;
304 default:
305 return getDkCardinality(operatorType, 3)*cardinality;
306 }
307 }
308
317 BasisPtr<DeviceType,outputValueType,pointValueType>
318 getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override{
319 if(subCellDim == 1) {
320 return Teuchos::rcp(new
322 (this->basisDegree_, pointType_));
323 } else if(subCellDim == 2) {
324 return Teuchos::rcp(new
326 (this->basisDegree_, pointType_));
327 }
328
329 INTREPID2_TEST_FOR_EXCEPTION(true,std::invalid_argument,"Input parameters out of bounds");
330 }
331
336 };
337
338}// namespace Intrepid2
339
341
342#endif
Header file for the abstract base class Intrepid2::Basis.
void getValues_HGRAD_Args(const outputValueViewType outputValues, const inputPointViewType inputPoints, const EOperator operatorType, const shards::CellTopology cellTopo, const ordinal_type basisCard)
Runtime check of the arguments for the getValues method in an HGRAD-conforming FEM basis....
Teuchos::RCP< Basis< DeviceType, OutputType, PointType > > BasisPtr
Basis Pointer.
Definition file for FEM basis functions of degree n for H(grad) functions on TET cells.
Header file for the Intrepid2::Basis_HGRAD_TET_Cn_FEM_ORTH class.
Header file for the Intrepid2::Basis_HGRAD_TRI_Cn_FEM class.
Header file for Intrepid2::PointTools class to provide utilities for barycentric coordinates,...
Implementation of the locally H(grad)-compatible FEM basis of variable order on the [-1,...
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Tetrahedron ce...
Kokkos::DynRankView< scalarType, DeviceType > vinv_
inverse of Generalized Vandermonde matrix, whose columns store the expansion coefficients of the noda...
EPointType pointType_
type of lattice used for creating the DoF coordinates
virtual const char * getName() const override
Returns basis name.
virtual bool requireOrientation() const override
True if orientation is required.
virtual void getDofCoeffs(ScalarViewType dofCoeffs) const override
Coefficients for computing degrees of freedom for Lagrangian basis If P is an element of the space sp...
BasisPtr< DeviceType, outputValueType, pointValueType > getSubCellRefBasis(const ordinal_type subCellDim, const ordinal_type subCellOrd) const override
returns the basis associated to a subCell.
virtual void getDofCoords(ScalarViewType dofCoords) const override
Returns spatial locations (coordinates) of degrees of freedom on the reference cell.
virtual void getValues(OutputViewType outputValues, const PointViewType inputPoints, const EOperator operatorType=OPERATOR_VALUE) const override
Evaluation of a FEM basis on a reference cell.
BasisPtr< typename Kokkos::HostSpace::device_type, outputValueType, pointValueType > getHostBasis() const override
Creates and returns a Basis object whose DeviceType template argument is Kokkos::HostSpace::device_ty...
Basis_HGRAD_TET_Cn_FEM(const ordinal_type order, const EPointType pointType=POINTTYPE_EQUISPACED)
Constructor.
Implementation of the default H(grad)-compatible Lagrange basis of arbitrary degree on Triangle cell.
An abstract base class that defines interface for concrete basis implementations for Finite Element (...
Kokkos::DynRankView< PointValueType, Kokkos::LayoutStride, DeviceType > PointViewType
View type for input points.
Kokkos::DynRankView< OutputValueType, Kokkos::LayoutStride, DeviceType > OutputViewType
View type for basis value output.
Kokkos::View< ordinal_type ***, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray3DHost
View type for 3d host array.
ordinal_type basisDegree_
Degree of the largest complete polynomial space that can be represented by the basis.
ordinal_type getCardinality() const
Returns cardinality of the basis.
Kokkos::DynRankView< scalarType, Kokkos::LayoutStride, DeviceType > ScalarViewType
View type for scalars.
Device DeviceType
(Kokkos) Device type on which Basis is templated. Does not necessarily return true for Kokkos::is_dev...
Kokkos::DynRankView< scalarType, DeviceType > dofCoords_
Coordinates of degrees-of-freedom for basis functions defined in physical space.
Kokkos::View< ordinal_type **, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray2DHost
View type for 2d host array.
ScalarTraits< pointValueType >::scalar_type scalarType
Scalar type for point values.
shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation https://trili...
Kokkos::View< ordinal_type *, typename ExecutionSpace::array_layout, Kokkos::HostSpace > OrdinalTypeArray1DHost
View type for 1d host array.
See Intrepid2::Basis_HGRAD_TET_Cn_FEM.
static constexpr ordinal_type MaxNumPtsPerBasisEval
The maximum number of points to eval in serial mode.
small utility functions