Intrepid2
Intrepid2_ProjectionToolsDefHCURL.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_PROJECTIONTOOLSDEFHCURL_HPP__
50 #define __INTREPID2_PROJECTIONTOOLSDEFHCURL_HPP__
51 
52 #include "Intrepid2_DefaultCubatureFactory.hpp"
53 #include "Intrepid2_ArrayTools.hpp"
54 #include "Intrepid2_FunctionSpaceTools.hpp"
55 
56 
57 namespace Intrepid2 {
58 namespace Experimental {
59 
60 template<typename SpT>
61 template<typename BasisType,
62 typename ortValueType, class ...ortProperties>
63 void
64 ProjectionTools<SpT>::getHCurlEvaluationPoints(typename BasisType::scalarViewType evaluationPoints,
65  typename BasisType::scalarViewType extDerivEvaluationPoints,
66  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
67  const BasisType* cellBasis,
68  ProjectionStruct<SpT, typename BasisType::scalarType> * projStruct,
69  const EvalPointsType evalPointType) {
70  typedef typename BasisType::scalarType scalarType;
71  typedef Kokkos::DynRankView<scalarType,SpT> scalarViewType;
72  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
73  const auto cellTopo = cellBasis->getBaseCellTopology();
74  ordinal_type dim = cellTopo.getDimension();
75  ordinal_type numCells = evaluationPoints.extent(0);
76  const ordinal_type edgeDim = 1;
77  const ordinal_type faceDim = 2;
78 
79  ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
80  ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
81 
82  Kokkos::View<ordinal_type*> eOrt("eOrt", numEdges), fOrt("fOrt", numFaces);
83 
84  for(ordinal_type ie=0; ie<numEdges; ++ie) {
85  range_type edgePointsRange;
86  scalarViewType cubPoints;
87  if(evalPointType == TARGET) {
88  edgePointsRange = projStruct->getTargetPointsRange(edgeDim, ie);
89  cubPoints = projStruct->getTargetEvalPoints(edgeDim, ie);
90  }
91  else {
92  edgePointsRange = projStruct->getBasisPointsRange(edgeDim, ie);
93  cubPoints = projStruct->getBasisEvalPoints(edgeDim, ie);
94  }
95 
96  scalarViewType orientedTargetCubPoints("orientedTargetCubPoints", cubPoints.extent(0),edgeDim);
97 
98  const auto topoKey = projStruct->getTopologyKey(edgeDim,ie);
99 
100  for(ordinal_type ic=0; ic<numCells; ++ic) {
101  orts(ic).getEdgeOrientation(eOrt.data(), numEdges);
102  ordinal_type ort = eOrt(ie);
103  Impl::OrientationTools::mapToModifiedReference(orientedTargetCubPoints,cubPoints,topoKey,ort);
104  CellTools<SpT>::mapToReferenceSubcell(Kokkos::subview(evaluationPoints,ic,edgePointsRange,Kokkos::ALL()), orientedTargetCubPoints, edgeDim, ie, cellBasis->getBaseCellTopology());
105  }
106  }
107 
108 
109  for(ordinal_type iface=0; iface<numFaces; ++iface) {
110 
111  scalarViewType cubPoints;//("cubPoints", numTargetCubPoints, faceDim);
112  range_type facePointsRange;
113  if(evalPointType == TARGET) {
114  cubPoints = projStruct->getTargetEvalPoints(faceDim, iface);
115  facePointsRange = projStruct->getTargetPointsRange(faceDim, iface);
116  } else {
117  cubPoints = projStruct->getBasisEvalPoints(faceDim, iface);
118  facePointsRange = projStruct->getBasisPointsRange(faceDim, iface);
119  }
120 
121  scalarViewType curlCubPoints;//("curlCubPoints", numTargetCurlCubPoints, faceDim);
122  range_type faceCurlPointsRange;
123  if(evalPointType == TARGET) {
124  curlCubPoints = projStruct->getTargetDerivEvalPoints(faceDim, iface);
125  faceCurlPointsRange = projStruct->getTargetDerivPointsRange(faceDim, iface);
126  } else {
127  curlCubPoints = projStruct->getBasisDerivEvalPoints(faceDim, iface);
128  faceCurlPointsRange = projStruct->getBasisDerivPointsRange(faceDim, iface);
129  }
130 
131  scalarViewType faceCubPoints("faceCubPoints", cubPoints.extent(0), faceDim);
132  scalarViewType faceCurlCubPoints("faceCurlCubPoints", curlCubPoints.extent(0), faceDim);
133 
134  const auto topoKey = projStruct->getTopologyKey(faceDim,iface);
135  for(ordinal_type ic=0; ic<numCells; ++ic) {
136 
137  orts(ic).getFaceOrientation(fOrt.data(), numFaces);
138 
139  ordinal_type ort = fOrt(iface);
140  Impl::OrientationTools::mapToModifiedReference(faceCubPoints,cubPoints,topoKey,ort);
141  CellTools<SpT>::mapToReferenceSubcell(Kokkos::subview(evaluationPoints, ic, facePointsRange, Kokkos::ALL()), faceCubPoints, faceDim, iface, cellBasis->getBaseCellTopology());
142 
143  Impl::OrientationTools::mapToModifiedReference(faceCurlCubPoints,curlCubPoints,topoKey,ort);
144  CellTools<SpT>::mapToReferenceSubcell(Kokkos::subview(extDerivEvaluationPoints, ic, faceCurlPointsRange, Kokkos::ALL()), faceCurlCubPoints, faceDim, iface, cellBasis->getBaseCellTopology());
145  }
146  }
147 
148 
149  if(cellBasis->getDofCount(dim,0)>0) {
150  range_type cellPointsRange;
151  scalarViewType cubPoints;
152  if(evalPointType == TARGET) {
153  cubPoints = projStruct->getTargetEvalPoints(dim, 0);
154  cellPointsRange = projStruct->getTargetPointsRange(dim, 0);
155  } else {
156  cubPoints = projStruct->getBasisEvalPoints(dim, 0);
157  cellPointsRange = projStruct->getBasisPointsRange(dim, 0);
158  }
159  RealSpaceTools<SpT>::clone(Kokkos::subview(evaluationPoints, Kokkos::ALL(), cellPointsRange, Kokkos::ALL()), cubPoints);
160 
161  range_type cellCurlPointsRange;
162  scalarViewType curlCubPoints;
163  if(evalPointType == TARGET) {
164  curlCubPoints = projStruct->getTargetDerivEvalPoints(dim, 0);
165  cellCurlPointsRange = projStruct->getTargetDerivPointsRange(dim, 0);
166  } else {
167  curlCubPoints = projStruct->getBasisDerivEvalPoints(dim, 0);
168  cellCurlPointsRange = projStruct->getBasisDerivPointsRange(dim, 0);
169  }
170  RealSpaceTools<SpT>::clone(Kokkos::subview(extDerivEvaluationPoints, Kokkos::ALL(), cellCurlPointsRange, Kokkos::ALL()), curlCubPoints);
171  }
172 }
173 
174 
175 template<typename SpT>
176 template<typename basisCoeffsValueType, class ...basisCoeffsProperties,
177 typename funValsValueType, class ...funValsProperties,
178 typename BasisType,
179 typename ortValueType,class ...ortProperties>
180 void
181 ProjectionTools<SpT>::getHCurlBasisCoeffs(Kokkos::DynRankView<basisCoeffsValueType,basisCoeffsProperties...> basisCoeffs,
182  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetAtEvalPoints,
183  const Kokkos::DynRankView<funValsValueType,funValsProperties...> targetCurlAtCurlEvalPoints,
184  const typename BasisType::scalarViewType evaluationPoints,
185  const typename BasisType::scalarViewType extDerivEvaluationPoints,
186  const Kokkos::DynRankView<ortValueType, ortProperties...> orts,
187  const BasisType* cellBasis,
188  ProjectionStruct<SpT, typename BasisType::scalarType> * projStruct){
189 
190  typedef typename Kokkos::Impl::is_space<SpT>::host_mirror_space::execution_space host_space_type;
191  typedef typename BasisType::scalarType scalarType;
192  typedef Kokkos::DynRankView<scalarType,SpT> scalarViewType;
193  typedef Kokkos::pair<ordinal_type,ordinal_type> range_type;
194  const auto cellTopo = cellBasis->getBaseCellTopology();
195  ordinal_type dim = cellTopo.getDimension();
196  ordinal_type numTotalEvaluationPoints(targetAtEvalPoints.extent(1)),
197  numTotalCurlEvaluationPoints(targetCurlAtCurlEvalPoints.extent(1));
198  ordinal_type basisCardinality = cellBasis->getCardinality();
199  ordinal_type numCells = targetAtEvalPoints.extent(0);
200  const ordinal_type edgeDim = 1;
201  const ordinal_type faceDim = 2;
202  const ordinal_type derDim = dim == 3 ? dim : 1;
203 
204  const std::string& name = cellBasis->getName();
205 
206  ordinal_type numEdges = (cellBasis->getDofCount(1, 0) > 0) ? cellTopo.getEdgeCount() : 0;
207  ordinal_type numFaces = (cellBasis->getDofCount(2, 0) > 0) ? cellTopo.getFaceCount() : 0;
208 
209  Kokkos::View<ordinal_type*> eOrt("eOrt", numEdges);
210  Kokkos::View<ordinal_type*> fOrt("fOrt", numFaces);
211  scalarViewType refEdgeTan("refEdgeTan", dim);
212  scalarViewType refFaceTangents("refFaceTangents", dim, 2);
213  scalarViewType refFaceNormal("refFaceNormal", dim);
214  auto refFaceTanU = Kokkos::subview(refFaceTangents, Kokkos::ALL, 0);
215  auto refFaceTanV = Kokkos::subview(refFaceTangents, Kokkos::ALL, 1);
216 
217  ordinal_type numEdgeDofs(0);
218  for(ordinal_type ie=0; ie<numEdges; ++ie)
219  numEdgeDofs += cellBasis->getDofCount(edgeDim,ie);
220 
221  ordinal_type numFaceDofs(0);
222  for(ordinal_type iface=0; iface<numFaces; ++iface)
223  numFaceDofs += cellBasis->getDofCount(faceDim,iface);
224 
225  Kokkos::View<ordinal_type*> computedDofs("computedDofs",numEdgeDofs+numFaceDofs);
226 
227  ordinal_type computedDofsCount = 0;
228 
229  ordinal_type numTotalCubPoints = projStruct->getNumBasisEvalPoints(), numTotalCurlCubPoints = projStruct->getNumBasisDerivEvalPoints();
230  scalarViewType cubPoints("cubPoints",numCells,numTotalCubPoints, dim);
231  scalarViewType curlCubPoints("curlCubPoints",numCells,numTotalCurlCubPoints, dim);
232  getHCurlEvaluationPoints(cubPoints, curlCubPoints, orts, cellBasis, projStruct, BASIS);
233 
234  scalarViewType basisAtCubPoints("basisAtCubPoints",numCells,basisCardinality, numTotalCubPoints, dim);
235  scalarViewType basisAtTargetCubPoints("basisAtTargetCubPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
236  {
237  scalarViewType nonOrientedBasisAtCubPoints("nonOrientedBasisAtCubPoints",numCells,basisCardinality, numTotalCubPoints, dim);
238  scalarViewType nonOrientedBasisAtTargetCubPoints("nonOrientedBasisAtTargetCubPoints",numCells,basisCardinality, numTotalEvaluationPoints, dim);
239  for(ordinal_type ic=0; ic<numCells; ++ic) {
240  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtTargetCubPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(evaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
241  cellBasis->getValues(Kokkos::subview(nonOrientedBasisAtCubPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(cubPoints, ic, Kokkos::ALL(), Kokkos::ALL()));
242  }
243 
244  OrientationTools<SpT>::modifyBasisByOrientation(basisAtCubPoints, nonOrientedBasisAtCubPoints, orts, cellBasis);
245  OrientationTools<SpT>::modifyBasisByOrientation(basisAtTargetCubPoints, nonOrientedBasisAtTargetCubPoints, orts, cellBasis);
246  }
247 
248  scalarViewType basisCurlAtCurlCubPoints;
249  scalarViewType basisCurlAtTargetCurlCubPoints;
250  if(numTotalCurlEvaluationPoints>0) {
251  scalarViewType nonOrientedBasisCurlAtTargetCurlCubPoints, nonOrientedBasisCurlAtCurlCubPoints;
252  if (dim == 3) {
253  basisCurlAtCurlCubPoints = scalarViewType ("basisCurlAtCurlCubPoints",numCells,basisCardinality, numTotalCurlCubPoints, dim);
254  nonOrientedBasisCurlAtCurlCubPoints = scalarViewType ("nonOrientedBasisCurlAtCurlCubPoints",numCells,basisCardinality, numTotalCurlCubPoints, dim);
255  basisCurlAtTargetCurlCubPoints = scalarViewType("basisCurlAtTargetCurlCubPoints",numCells,basisCardinality, numTotalCurlEvaluationPoints, dim);
256  nonOrientedBasisCurlAtTargetCurlCubPoints = scalarViewType("nonOrientedBasisCurlAtTargetCurlCubPoints",numCells,basisCardinality, numTotalCurlEvaluationPoints, dim);
257  } else {
258  basisCurlAtCurlCubPoints = scalarViewType ("basisCurlAtCurlCubPoints",numCells,basisCardinality, numTotalCurlCubPoints);
259  nonOrientedBasisCurlAtCurlCubPoints = scalarViewType ("nonOrientedBasisCurlAtCurlCubPoints",numCells,basisCardinality, numTotalCurlCubPoints);
260  basisCurlAtTargetCurlCubPoints = scalarViewType("basisCurlAtTargetCurlCubPoints",numCells,basisCardinality, numTotalCurlEvaluationPoints);
261  nonOrientedBasisCurlAtTargetCurlCubPoints = scalarViewType("nonOrientedBasisCurlAtTargetCurlCubPoints",numCells,basisCardinality, numTotalCurlEvaluationPoints);
262  }
263  for(ordinal_type ic=0; ic<numCells; ++ic) {
264  cellBasis->getValues(Kokkos::subview(nonOrientedBasisCurlAtCurlCubPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(curlCubPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_CURL);
265  cellBasis->getValues(Kokkos::subview(nonOrientedBasisCurlAtTargetCurlCubPoints,ic,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()), Kokkos::subview(extDerivEvaluationPoints, ic, Kokkos::ALL(), Kokkos::ALL()),OPERATOR_CURL);
266  }
267  OrientationTools<SpT>::modifyBasisByOrientation(basisCurlAtCurlCubPoints, nonOrientedBasisCurlAtCurlCubPoints, orts, cellBasis);
268  OrientationTools<SpT>::modifyBasisByOrientation(basisCurlAtTargetCurlCubPoints, nonOrientedBasisCurlAtTargetCurlCubPoints, orts, cellBasis);
269  }
270 
271  for(ordinal_type ie=0; ie<numEdges; ++ie) {
272 
273  ordinal_type edgeCardinality = cellBasis->getDofCount(edgeDim,ie);
274  ordinal_type numCubPoints = projStruct->getNumBasisEvalPoints(edgeDim, ie);
275  ordinal_type numTargetCubPoints = projStruct->getNumTargetEvalPoints(edgeDim, ie);
276 
277  CellTools<SpT>::getReferenceEdgeTangent(refEdgeTan, ie, cellBasis->getBaseCellTopology());
278 
279  scalarViewType tanBasisAtElemCubPoints("tanBasisAtElemCubPoints",numCells,edgeCardinality, numCubPoints);
280  scalarViewType tanBasisAtTargetCubPoints("tanBasisAtTargetCubPoints",numCells,edgeCardinality, numTargetCubPoints);
281  scalarViewType weightedTanBasisAtElemCubPoints("weightedTanBasisAtElemCubPoints",numCells,edgeCardinality, numCubPoints);
282  scalarViewType weightedTanBasisAtTargetCubPoints("weightedTanBasisAtTargetCubPoints",numCells,edgeCardinality, numTargetCubPoints);
283  scalarViewType tanTargetAtTargetCubPoints("normalTargetAtTargetCubPoints",numCells, numTargetCubPoints);
284 
285  scalarViewType targetEvalWeights = projStruct->getTargetEvalWeights(edgeDim, ie);
286  scalarViewType basisEvalWeights = projStruct->getBasisEvalWeights(edgeDim, ie);
287 
288  //Note: we are not considering the jacobian of the orientation map since it is simply a scalar term for the integrals and it does not affect the projection
289  ordinal_type offsetBasis = projStruct->getBasisPointsRange(edgeDim, ie).first;
290  ordinal_type offsetTarget = projStruct->getTargetPointsRange(edgeDim, ie).first;
291  for(ordinal_type ic=0; ic<numCells; ++ic) {
292  for(ordinal_type j=0; j <edgeCardinality; ++j) {
293  ordinal_type jdof = cellBasis->getDofOrdinal(edgeDim, ie, j);
294  for(ordinal_type iq=0; iq <numCubPoints; ++iq) {
295  for(ordinal_type d=0; d <dim; ++d)
296  tanBasisAtElemCubPoints(ic,j,iq) += refEdgeTan(d)*basisAtCubPoints(ic,jdof,offsetBasis+iq,d);
297  weightedTanBasisAtElemCubPoints(ic,j,iq) = tanBasisAtElemCubPoints(ic,j,iq)*basisEvalWeights(iq);
298  }
299  for(ordinal_type iq=0; iq <numTargetCubPoints; ++iq) {
300  for(ordinal_type d=0; d <dim; ++d)
301  tanBasisAtTargetCubPoints(ic,j,iq) += refEdgeTan(d)*basisAtTargetCubPoints(ic,jdof,offsetTarget+iq,d);
302  weightedTanBasisAtTargetCubPoints(ic,j,iq) = tanBasisAtTargetCubPoints(ic,j,iq)*targetEvalWeights(iq);
303  }
304  }
305  for(ordinal_type iq=0; iq <numTargetCubPoints; ++iq)
306  for(ordinal_type d=0; d <dim; ++d)
307  tanTargetAtTargetCubPoints(ic,iq) += refEdgeTan(d)*targetAtEvalPoints(ic,offsetTarget+iq,d);
308  }
309  scalarViewType edgeMassMat_("edgeMassMat_", numCells, edgeCardinality+1, edgeCardinality+1),
310  edgeRhsMat_("rhsMat_", numCells, edgeCardinality+1);
311 
312  scalarViewType cubWeights_("cubWeights_", numCells, 1, basisEvalWeights.extent(0)), targetEvalWeights_("targetEvalWeights", numCells, 1, targetEvalWeights.extent(0));
313  RealSpaceTools<SpT>::clone(cubWeights_, basisEvalWeights);
314  RealSpaceTools<SpT>::clone(targetEvalWeights_, targetEvalWeights);
315 
316  range_type range_H(0, edgeCardinality);
317  range_type range_B(edgeCardinality, edgeCardinality+1);
318  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(edgeMassMat_,Kokkos::ALL(),range_H,range_H), tanBasisAtElemCubPoints, weightedTanBasisAtElemCubPoints);
319  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(edgeMassMat_,Kokkos::ALL(),range_H,range_B), tanBasisAtElemCubPoints, cubWeights_);
320  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(edgeRhsMat_,Kokkos::ALL(),range_H), tanTargetAtTargetCubPoints, weightedTanBasisAtTargetCubPoints);
321  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(edgeRhsMat_,Kokkos::ALL(),range_B), tanTargetAtTargetCubPoints, targetEvalWeights_);
322  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> edgeMassMat("edgeMassMat", edgeCardinality+1,edgeCardinality+1);
323  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> edgeRhsMat("edgeRhsMat",edgeCardinality+1, 1);
324 
325  Teuchos::LAPACK<ordinal_type,funValsValueType> lapack;
326  ordinal_type info = 0;
327  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> pivVec("pivVec", edgeCardinality+1, 1);
328  for(ordinal_type ic=0; ic<numCells; ++ic) {
329  Kokkos::deep_copy(edgeMassMat,funValsValueType(0)); //LAPACK might overwrite the matrix
330 
331  for(ordinal_type i=0; i<edgeCardinality; ++i) {
332  edgeRhsMat(i,0) = edgeRhsMat_(ic,i);
333  for(ordinal_type j=0; j<edgeCardinality+1; ++j)
334  edgeMassMat(i,j) = edgeMassMat_(ic,i,j);
335  edgeMassMat(edgeCardinality,i) = edgeMassMat_(ic,i,edgeCardinality);
336  }
337  edgeRhsMat(edgeCardinality,0) = edgeRhsMat_(ic,edgeCardinality);
338 
339  lapack.GESV(edgeCardinality+1, 1,
340  edgeMassMat.data(),
341  edgeMassMat.stride_1(),
342  (ordinal_type*)pivVec.data(),
343  edgeRhsMat.data(),
344  edgeRhsMat.stride_1(),
345  &info);
346 
347  if (info) {
348  std::stringstream ss;
349  ss << ">>> ERROR (Intrepid::ProjectionTools::getBasisCoeffs): "
350  << "LAPACK return with error code: "
351  << info;
352  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
353  }
354 
355  for(ordinal_type i=0; i<edgeCardinality; ++i){
356  ordinal_type edge_dof = cellBasis->getDofOrdinal(edgeDim, ie, i);
357  basisCoeffs(ic,edge_dof) = edgeRhsMat(i,0);
358  }
359  }
360  for(ordinal_type i=0; i<edgeCardinality; ++i)
361  computedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(edgeDim, ie, i);
362 
363  }
364 
365  scalarViewType ortJacobian("ortJacobian", faceDim, faceDim);
366 
367  Basis<host_space_type,scalarType,scalarType> *hgradBasis = NULL;
368  for(ordinal_type iface=0; iface<numFaces; ++iface) {
369 
370  if(name.find("HEX")!=std::string::npos)
371  hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM<host_space_type,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
372  else if(name.find("TET")!=std::string::npos)
373  hgradBasis = new Basis_HGRAD_TRI_Cn_FEM<host_space_type,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
374  else {
375  std::stringstream ss;
376  ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): "
377  << "Method not implemented for basis " << name;
378  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
379  return;
380  }
381 
382 
383  ordinal_type numFaceDofs = cellBasis->getDofCount(faceDim,iface);
384  ordinal_type numTargetCubPoints = projStruct->getNumTargetEvalPoints(faceDim, iface);
385 
386  ordinal_type numTargetCurlCubPoints = projStruct->getNumTargetDerivEvalPoints(faceDim, iface);
387  ordinal_type numCubPoints = projStruct->getNumBasisEvalPoints(faceDim, iface);
388 
389  scalarViewType hgradBasisGradAtCubPoints("hgradBasisGradAtCubPoints",hgradBasis->getCardinality(), numCubPoints, faceDim);
390  scalarViewType hgradBasisGradAtTargetCubPoints("hgradBasisGradAtTargetCubPoints",hgradBasis->getCardinality(), numTargetCubPoints, faceDim);
391 
392  ordinal_type internalHgradCardinality = hgradBasis->getDofCount(faceDim,0);
393  scalarViewType internalHgradBasisGradAtCubPoints("internalHgradBasisGradAtCubPoints",1, internalHgradCardinality, numCubPoints, faceDim);
394  scalarViewType internalHgradBasisGradAtTargetCubPoints("internalHgradBasisGradAtTargetCubPoints",1, internalHgradCardinality, numTargetCubPoints, faceDim);
395 
396 
397  CellTools<SpT>::getReferenceFaceNormal(refFaceNormal, iface, cellTopo);
398  CellTools<SpT>::getReferenceFaceTangents(refFaceTanU, refFaceTanV,iface, cellTopo);
399 
400  hgradBasis->getValues(hgradBasisGradAtCubPoints,projStruct->getBasisEvalPoints(faceDim, iface), OPERATOR_GRAD);
401  hgradBasis->getValues(hgradBasisGradAtTargetCubPoints,projStruct->getTargetEvalPoints(faceDim, iface),OPERATOR_GRAD);
402 
403  for(ordinal_type j=0; j <internalHgradCardinality; ++j) {
404  ordinal_type face_dof = hgradBasis->getDofOrdinal(faceDim, 0, j);
405  for(ordinal_type d=0; d <faceDim; ++d) {
406  for(ordinal_type iq=0; iq <numCubPoints; ++iq)
407  internalHgradBasisGradAtCubPoints(0,j,iq,d)= hgradBasisGradAtCubPoints(face_dof,iq,d);
408  for(ordinal_type iq=0; iq <numTargetCubPoints; ++iq)
409  internalHgradBasisGradAtTargetCubPoints(0,j,iq,d)= hgradBasisGradAtTargetCubPoints(face_dof,iq,d);
410  }
411  }
412 
413  scalarViewType tanBasisAtElemCubPoints("tanBasisAtElemCubPoints",numCells,numFaceDofs, numCubPoints,dim-1);
414  scalarViewType tanBasisAtTargetCubPoints("tanBasisAtTargetCubPoints",numCells,numFaceDofs, numTargetCubPoints,dim-1);
415  scalarViewType normalBasisCurlAtElemCubPoints("normaBasisCurlAtElemCubPoints",numCells,numFaceDofs, numCubPoints);
416  scalarViewType wNormalBasisCurlAtElemCubPoints("weightedNormalBasisCurlAtElemCubPoints",numCells,numFaceDofs, numCubPoints);
417 
418  scalarViewType tanTargetAtTargetCubPoints("tanTargetAtTargetCubPoints",numCells, numTargetCubPoints, dim-1);
419  scalarViewType normalTargetCurlAtTargetCubPoints("normalTargetCurlAtTargetCubPoints",numCells, numTargetCurlCubPoints);
420  scalarViewType normalBasisCurlAtTargetCurlCubPoints("normalBasisCurlAtTargetCurlCubPoints",numCells,numFaceDofs, numTargetCurlCubPoints);
421  scalarViewType wNormalBasisCurlBasisAtTargetCurlCubPoints("weightedNormalBasisCurlAtTargetCurlCubPoints",numCells,numFaceDofs, numTargetCurlCubPoints);
422 
423  scalarViewType wHgradBasisGradAtCubPoints("wHgradBasisGradAtCubPoints",1, internalHgradCardinality, numCubPoints, faceDim);
424  scalarViewType wHgradBasisGradAtCubPoints_("wHgradBasisGradAtCubPoints_",numCells, internalHgradCardinality, numCubPoints, faceDim);
425  scalarViewType wHgradBasisGradAtTargetCubPoints("wHgradBasisGradAtTargetCubPoints",1, internalHgradCardinality, numTargetCubPoints, faceDim);
426  scalarViewType wHgradBasisGradAtTargetCubPoints_("wHgradBasisGradAtTargetCubPoints_",numCells, internalHgradCardinality, numTargetCubPoints, faceDim);
427 
428  scalarViewType mNormalComputedProjectionCurl("mNormalComputedProjection", numCells,numCubPoints);
429  scalarViewType mTanComputedProjection("mTanComputedProjection", numCells,numCubPoints,dim-1);
430 
431  scalarViewType targetDerivEvalWeights = projStruct->getTargetDerivEvalWeights(faceDim, iface);
432  ordinal_type offsetBasis = projStruct->getBasisPointsRange(faceDim, iface).first;
433  ordinal_type offsetBasisCurl = projStruct->getBasisDerivPointsRange(faceDim, iface).first;
434  ordinal_type offsetTarget = projStruct->getTargetPointsRange(faceDim, iface).first;
435  ordinal_type offsetTargetCurl = projStruct->getTargetDerivPointsRange(faceDim, iface).first;
436 
437 
438  //Note: we are not considering the jacobian of the orientation map since it is simply a scalar term for the integrals and it does not affect the projection
439  const auto topoKey = projStruct->getTopologyKey(faceDim,iface);
440  for(ordinal_type ic=0; ic<numCells; ++ic) {
441  orts(ic).getFaceOrientation(fOrt.data(), numFaces);
442 
443  ordinal_type ort = fOrt(iface);
444 
445  Impl::OrientationTools::getJacobianOfOrientationMap(ortJacobian, topoKey, ort);
446 
447  for(ordinal_type j=0; j <numFaceDofs; ++j) {
448  ordinal_type jdof = cellBasis->getDofOrdinal(faceDim, iface, j);
449  for(ordinal_type iq=0; iq <numCubPoints; ++iq)
450  for(ordinal_type d=0; d <dim; ++d) {
451  normalBasisCurlAtElemCubPoints(ic,j,iq) += refFaceNormal(d)*basisCurlAtCurlCubPoints(ic,jdof,offsetBasisCurl+iq,d);
452  for(ordinal_type itan=0; itan <dim-1; ++itan)
453  for(ordinal_type jtan=0; jtan <dim-1; ++jtan)
454  tanBasisAtElemCubPoints(ic,j,iq,itan) += refFaceTangents(d,jtan)*ortJacobian(jtan,itan)*basisAtCubPoints(ic,jdof,offsetBasis+iq,d);
455  }
456  for(ordinal_type iq=0; iq <numTargetCurlCubPoints; ++iq) {
457  for(ordinal_type d=0; d <dim; ++d)
458  normalBasisCurlAtTargetCurlCubPoints(ic,j,iq) += refFaceNormal(d)*basisCurlAtTargetCurlCubPoints(ic,jdof,offsetTargetCurl+iq,d);
459  wNormalBasisCurlBasisAtTargetCurlCubPoints(ic,j,iq) = normalBasisCurlAtTargetCurlCubPoints(ic,j,iq)*targetDerivEvalWeights[iq];
460  }
461  }
462  for(ordinal_type j=0; j <numEdgeDofs; ++j) {
463  ordinal_type jdof = computedDofs(j);
464  for(ordinal_type iq=0; iq <numCubPoints; ++iq)
465  for(ordinal_type d=0; d <dim; ++d) {
466  mNormalComputedProjectionCurl(ic,iq) -= refFaceNormal(d)*basisCoeffs(ic,jdof)*basisCurlAtCurlCubPoints(ic,jdof,offsetBasisCurl+iq,d);
467  for(ordinal_type itan=0; itan <dim-1; ++itan)
468  for(ordinal_type jtan=0; jtan <dim-1; ++jtan)
469  mTanComputedProjection(ic,iq,itan) -= refFaceTangents(d,jtan)*ortJacobian(jtan,itan)*basisCoeffs(ic,jdof)*basisAtCubPoints(ic,jdof,offsetBasis+iq,d);
470  }
471  }
472  for(ordinal_type iq=0; iq <numTargetCubPoints; ++iq)
473  for(ordinal_type d=0; d <dim; ++d)
474  for(ordinal_type itan=0; itan <dim-1; ++itan)
475  for(ordinal_type jtan=0; jtan <dim-1; ++jtan)
476  tanTargetAtTargetCubPoints(ic,iq,itan) += refFaceTangents(d,jtan)*ortJacobian(jtan,itan)*targetAtEvalPoints(ic,offsetTarget+iq,d);
477  for(ordinal_type iq=0; iq <numTargetCurlCubPoints; ++iq)
478  for(ordinal_type d=0; d <dim; ++d)
479  normalTargetCurlAtTargetCubPoints(ic,iq) += refFaceNormal(d)*targetCurlAtCurlEvalPoints(ic,offsetTargetCurl+iq,d);
480  }
481 
482  scalarViewType faceMassMat_("faceMassMat_", numCells, numFaceDofs+internalHgradCardinality, numFaceDofs+internalHgradCardinality),
483  faceRhsMat_("rhsMat_", numCells, numFaceDofs+internalHgradCardinality);
484 
485  scalarViewType targetCubWeights_("targetCubWeights_", 1, projStruct->getNumTargetEvalPoints(faceDim, iface));
486  RealSpaceTools<SpT>::clone(targetCubWeights_, projStruct->getTargetEvalWeights(faceDim, iface));
487  scalarViewType cubWeights_("cubWeights_", numCells, 1, numCubPoints);
488  RealSpaceTools<SpT>::clone(cubWeights_, projStruct->getBasisEvalWeights(faceDim, iface));
489  ArrayTools<SpT>::scalarMultiplyDataField( wNormalBasisCurlAtElemCubPoints, Kokkos::subview(cubWeights_, Kokkos::ALL(),0, Kokkos::ALL()),normalBasisCurlAtElemCubPoints, false);
490  ArrayTools<SpT>::scalarMultiplyDataField( wHgradBasisGradAtCubPoints, Kokkos::subview(cubWeights_, 0, Kokkos::ALL(), Kokkos::ALL()),internalHgradBasisGradAtCubPoints, false);
491  ArrayTools<SpT>::scalarMultiplyDataField( wHgradBasisGradAtTargetCubPoints, targetCubWeights_, internalHgradBasisGradAtTargetCubPoints , false);
492 
493  RealSpaceTools<SpT>::clone(wHgradBasisGradAtCubPoints_,Kokkos::subview(wHgradBasisGradAtCubPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
494  RealSpaceTools<SpT>::clone(wHgradBasisGradAtTargetCubPoints_,Kokkos::subview(wHgradBasisGradAtTargetCubPoints,0,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL()));
495 
496  range_type range_H(0, numFaceDofs);
497  range_type range_B(numFaceDofs, numFaceDofs+internalHgradCardinality);
498  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_H), normalBasisCurlAtElemCubPoints, wNormalBasisCurlAtElemCubPoints);
499  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(faceMassMat_,Kokkos::ALL(),range_H,range_B), tanBasisAtElemCubPoints, wHgradBasisGradAtCubPoints_);
500 
501  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), normalTargetCurlAtTargetCubPoints, wNormalBasisCurlBasisAtTargetCurlCubPoints);
502  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_H), mNormalComputedProjectionCurl, wNormalBasisCurlAtElemCubPoints,true);
503 
504  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_B), tanTargetAtTargetCubPoints, wHgradBasisGradAtTargetCubPoints_);
505  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(faceRhsMat_,Kokkos::ALL(),range_B), mTanComputedProjection, wHgradBasisGradAtCubPoints_,true);
506 
507  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> faceMassMat("faceMassMat", numFaceDofs+internalHgradCardinality,numFaceDofs+internalHgradCardinality);
508  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> faceRhsMat("faceRhsMat",numFaceDofs+internalHgradCardinality, 1);
509 
510  Teuchos::LAPACK<ordinal_type,funValsValueType> lapack;
511  ordinal_type info = 0;
512  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> pivVec("pivVec", numFaceDofs+internalHgradCardinality, 1);
513  for(ordinal_type ic=0; ic<numCells; ++ic) {
514  Kokkos::deep_copy(faceMassMat,funValsValueType(0)); //LAPACK might overwrite the matrix
515 
516  for(ordinal_type i=0; i<numFaceDofs; ++i) {
517  for(ordinal_type j=0; j<numFaceDofs+internalHgradCardinality; ++j)
518  faceMassMat(i,j) = faceMassMat_(ic,i,j);
519  for(ordinal_type j=0; j<internalHgradCardinality; ++j)
520  faceMassMat(numFaceDofs+j,i) = faceMassMat_(ic,i,numFaceDofs+j);
521  }
522 
523  for(ordinal_type i=0; i<numFaceDofs+internalHgradCardinality; ++i)
524  faceRhsMat(i,0) = faceRhsMat_(ic,i);
525 
526  lapack.GESV(numFaceDofs+internalHgradCardinality, 1,
527  faceMassMat.data(),
528  faceMassMat.stride_1(),
529  (ordinal_type*)pivVec.data(),
530  faceRhsMat.data(),
531  faceRhsMat.stride_1(),
532  &info);
533 
534  if (info) {
535  std::stringstream ss;
536  ss << ">>> ERROR (Intrepid::ProjectionTools::getBasisCoeffs): "
537  << "LAPACK return with error code: "
538  << info;
539  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
540  }
541 
542  for(ordinal_type i=0; i<numFaceDofs; ++i){
543  ordinal_type face_dof = cellBasis->getDofOrdinal(faceDim, iface, i);
544  basisCoeffs(ic,face_dof) = faceRhsMat(i,0);
545  }
546  }
547 
548  for(ordinal_type i=0; i<numFaceDofs; ++i)
549  computedDofs(computedDofsCount++) = cellBasis->getDofOrdinal(faceDim, iface, i);
550 
551  delete hgradBasis;
552  }
553 
554  ordinal_type numElemDofs = cellBasis->getDofCount(dim,0);
555  if(numElemDofs>0) {
556  if(name.find("HEX")!=std::string::npos)
557  hgradBasis = new Basis_HGRAD_HEX_Cn_FEM<host_space_type,scalarType,scalarType>(cellBasis->getDegree());
558  else if(name.find("TET")!=std::string::npos)
559  hgradBasis = new Basis_HGRAD_TET_Cn_FEM<host_space_type,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
560  else if(name.find("TRI")!=std::string::npos)
561  hgradBasis = new Basis_HGRAD_TRI_Cn_FEM<host_space_type,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
562  else if(name.find("QUAD")!=std::string::npos)
563  hgradBasis = new Basis_HGRAD_QUAD_Cn_FEM<host_space_type,scalarType,scalarType>(cellBasis->getDegree(),POINTTYPE_WARPBLEND);
564  else {
565  std::stringstream ss;
566  ss << ">>> ERROR (Intrepid2::ProjectionTools::getHCurlBasisCoeffs): "
567  << "Method not implemented for basis " << name;
568  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
569  return;
570  }
571 
572  range_type cellPointsRange = projStruct->getTargetPointsRange(dim, 0);
573  range_type cellCurlPointsRange = projStruct->getTargetDerivPointsRange(dim, 0);
574 
575  ordinal_type numTargetCurlCubPoints = projStruct->getNumTargetDerivEvalPoints(dim,0);
576  ordinal_type numCubPoints = projStruct->getNumBasisEvalPoints(dim,0);
577  ordinal_type numTargetCubPoints = projStruct->getNumTargetEvalPoints(dim,0);
578 
579  scalarViewType hgradBasisGradAtCubPoints("hgradBasisGradAtCubPoints",hgradBasis->getCardinality(), numCubPoints, dim);
580  scalarViewType hgradBasisGradAtTargetCubPoints("hgradBasisGradAtTargetCubPoints",hgradBasis->getCardinality(), numTargetCubPoints, dim);
581 
582  ordinal_type internalHgradCardinality = hgradBasis->getDofCount(dim,0);
583  scalarViewType internalHgradBasisGradAtCubPoints("internalHgradBasisGradAtCubPoints",1, internalHgradCardinality, numCubPoints, dim);
584  scalarViewType internalHgradBasisGradAtTargetCubPoints("internalHgradBasisGradAtTargetCubPoints",numCells, internalHgradCardinality, numTargetCubPoints, dim);
585  scalarViewType wHgradBasisGradAtTargetCubPoints("wHgradBasisGradAtTargetCubPoints",numCells, internalHgradCardinality, numTargetCubPoints, dim);
586  scalarViewType wHgradBasisGradAtCubPoints("wHgradBasisGradAtCubPoints",numCells, internalHgradCardinality, numCubPoints, dim);
587 
588  scalarViewType targetEvalWeights = projStruct->getTargetEvalWeights(dim, 0);
589  scalarViewType basisEvalWeights = projStruct->getBasisEvalWeights(dim, 0);
590 
591  hgradBasis->getValues(hgradBasisGradAtCubPoints,projStruct->getBasisEvalPoints(dim, 0), OPERATOR_GRAD);
592  hgradBasis->getValues(hgradBasisGradAtTargetCubPoints,projStruct->getTargetEvalPoints(dim, 0),OPERATOR_GRAD);
593 
594  for(ordinal_type j=0; j <internalHgradCardinality; ++j) {
595  ordinal_type idof = hgradBasis->getDofOrdinal(dim, 0, j);
596  for(ordinal_type d=0; d <dim; ++d) {
597  for(ordinal_type iq=0; iq <numCubPoints; ++iq) {
598  internalHgradBasisGradAtCubPoints(0,j,iq,d)= hgradBasisGradAtCubPoints(idof,iq,d);
599  for(ordinal_type ic=0; ic<numCells; ++ic)
600  wHgradBasisGradAtCubPoints(ic,j,iq,d) = internalHgradBasisGradAtCubPoints(0,j,iq,d)*basisEvalWeights[iq];
601  }
602  for(ordinal_type iq=0; iq <numTargetCubPoints; ++iq)
603  for(ordinal_type ic=0; ic<numCells; ++ic) {
604  internalHgradBasisGradAtTargetCubPoints(ic,j,iq,d)= hgradBasisGradAtTargetCubPoints(idof,iq,d);
605  wHgradBasisGradAtTargetCubPoints(ic,j,iq,d) = internalHgradBasisGradAtTargetCubPoints(ic,j,iq,d)*targetEvalWeights[iq];
606  }
607  }
608  }
609 
610  scalarViewType internalBasisAtElemcubPoints("basisElemAtcubPoints",numCells,numElemDofs, numCubPoints, dim);
611  scalarViewType internalBasisCurlAtElemcubPoints("internalBasisCurlAtElemcubPoints",numCells,numElemDofs, numCubPoints, derDim);
612  scalarViewType internalBasisCurlAtTargetCurlCubPoints("weightedBasisCurlAtElemCubPoints",numCells,numElemDofs, numTargetCurlCubPoints,derDim);
613  scalarViewType mComputedProjectionCurl("mComputedProjectionCurl", numCells, numCubPoints, derDim);
614  scalarViewType mComputedProjection("mComputedProjection", numCells, numCubPoints, dim);
615  scalarViewType wBasisCurlAtElemCubPoints("weightedBasisCurlAtElemCubPoints",numCells,numElemDofs, numCubPoints,derDim);
616  scalarViewType wBasisCurlBasisAtTargetCurlCubPoints("weightedBasisCurlAtTargetCurlCubPoints",numCells,numElemDofs, numTargetCurlCubPoints,derDim);
617  scalarViewType targetDerivEvalWeights = projStruct->getTargetDerivEvalWeights(dim, 0);
618  ordinal_type offsetBasis = projStruct->getBasisPointsRange(dim, 0).first;
619  ordinal_type offsetBasisCurl = projStruct->getBasisDerivPointsRange(dim, 0).first;
620  ordinal_type offsetTargetCurl = projStruct->getTargetDerivPointsRange(dim, 0).first;
621 
622  for(ordinal_type j=0; j <numElemDofs; ++j) {
623  ordinal_type idof = cellBasis->getDofOrdinal(dim, 0, j);
624  for(ordinal_type ic=0; ic<numCells; ++ic) {
625  for(ordinal_type d=0; d <dim; ++d)
626  for(ordinal_type iq=0; iq <numCubPoints; ++iq)
627  internalBasisAtElemcubPoints(ic,j,iq,d)=basisAtCubPoints(ic,idof,offsetBasis+iq,d);
628 
629  for(ordinal_type d=0; d <derDim; ++d) {
630  for(ordinal_type iq=0; iq <numCubPoints; ++iq) {
631  internalBasisCurlAtElemcubPoints(ic,j,iq,d)=basisCurlAtCurlCubPoints(ic,idof,offsetBasisCurl+iq,d);
632  wBasisCurlAtElemCubPoints(ic,j,iq,d)=internalBasisCurlAtElemcubPoints(ic,j,iq,d)*basisEvalWeights[iq];
633  }
634  for(ordinal_type iq=0; iq <numTargetCurlCubPoints; ++iq) {
635  internalBasisCurlAtTargetCurlCubPoints(ic,j,iq,d)= basisCurlAtTargetCurlCubPoints(ic,idof,offsetTargetCurl+iq,d);
636  wBasisCurlBasisAtTargetCurlCubPoints(ic,j,iq,d) = internalBasisCurlAtTargetCurlCubPoints(ic,j,iq,d)*targetDerivEvalWeights[iq];
637  }
638  }
639  }
640  }
641 
642  for(ordinal_type j=0; j <numEdgeDofs+numFaceDofs; ++j) {
643  ordinal_type jdof = computedDofs(j);
644  for(ordinal_type iq=0; iq <numCubPoints; ++iq)
645  for(ordinal_type ic=0; ic<numCells; ++ic) {
646  for(ordinal_type d=0; d <derDim; ++d)
647  mComputedProjectionCurl(ic,iq,d) -= basisCoeffs(ic,jdof)*basisCurlAtCurlCubPoints(ic,jdof,offsetBasisCurl+iq,d);
648  for(ordinal_type d=0; d <dim; ++d)
649  mComputedProjection(ic,iq,d) -= basisCoeffs(ic,jdof)*basisAtCubPoints(ic,jdof,offsetBasis+iq,d);
650  }
651  }
652 
653  scalarViewType cellMassMat_("cellMassMat_", numCells, numElemDofs+internalHgradCardinality, numElemDofs+internalHgradCardinality),
654  cellRhsMat_("rhsMat_", numCells, numElemDofs+internalHgradCardinality);
655 
656  range_type range_H(0, numElemDofs);
657  range_type range_B(numElemDofs, numElemDofs+internalHgradCardinality);
658  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_H), internalBasisCurlAtElemcubPoints, wBasisCurlAtElemCubPoints);
659  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellMassMat_,Kokkos::ALL(),range_H,range_B), internalBasisAtElemcubPoints, wHgradBasisGradAtCubPoints);
660  if(dim==3)
661  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), Kokkos::subview(targetCurlAtCurlEvalPoints,Kokkos::ALL(),cellCurlPointsRange,Kokkos::ALL()), wBasisCurlBasisAtTargetCurlCubPoints);
662  else
663  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), Kokkos::subview(targetCurlAtCurlEvalPoints,Kokkos::ALL(),cellCurlPointsRange), Kokkos::subview(wBasisCurlBasisAtTargetCurlCubPoints,Kokkos::ALL(),Kokkos::ALL(),Kokkos::ALL(),0));
664  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_H), mComputedProjectionCurl, wBasisCurlAtElemCubPoints, true);
665 
666  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_B), Kokkos::subview(targetAtEvalPoints,Kokkos::ALL(),cellPointsRange,Kokkos::ALL()), wHgradBasisGradAtTargetCubPoints);
667  FunctionSpaceTools<SpT >::integrate(Kokkos::subview(cellRhsMat_,Kokkos::ALL(),range_B), mComputedProjection, wHgradBasisGradAtCubPoints, true);
668 
669  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> cellMassMat("cellMassMat", numElemDofs+internalHgradCardinality,numElemDofs+internalHgradCardinality);
670  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> cellRhsMat("cellRhsMat",numElemDofs+internalHgradCardinality, 1);
671 
672  Teuchos::LAPACK<ordinal_type,funValsValueType> lapack;
673  ordinal_type info = 0;
674  Kokkos::View<funValsValueType**,Kokkos::LayoutLeft,host_space_type> pivVec("pivVec", numElemDofs+internalHgradCardinality, 1);
675 
676  for(ordinal_type ic=0; ic<numCells; ++ic) {
677  Kokkos::deep_copy(cellMassMat,funValsValueType(0)); //LAPACK might overwrite the matrix
678 
679  for(ordinal_type i=0; i<numElemDofs; ++i) {
680  for(ordinal_type j=0; j<numElemDofs+internalHgradCardinality; ++j)
681  cellMassMat(i,j) = cellMassMat_(ic,i,j);
682  for(ordinal_type j=0; j<internalHgradCardinality; ++j)
683  cellMassMat(numElemDofs+j,i) = cellMassMat_(0,i,numElemDofs+j);
684  }
685 
686  for(ordinal_type i=0; i<numElemDofs+internalHgradCardinality; ++i)
687  cellRhsMat(i,0) = cellRhsMat_(ic,i);
688 
689 
690  lapack.GESV(numElemDofs+internalHgradCardinality, 1,
691  cellMassMat.data(),
692  cellMassMat.stride_1(),
693  (ordinal_type*)pivVec.data(),
694  cellRhsMat.data(),
695  cellRhsMat.stride_1(),
696  &info);
697 
698  for(ordinal_type i=0; i<numElemDofs; ++i){
699  ordinal_type idof = cellBasis->getDofOrdinal(dim, 0, i);
700  basisCoeffs(ic,idof) = cellRhsMat(i,0);
701  }
702 
703  if (info) {
704  std::stringstream ss;
705  ss << ">>> ERROR (Intrepid::ProjectionTools::getBasisCoeffs): "
706  << "LAPACK return with error code: "
707  << info;
708  INTREPID2_TEST_FOR_EXCEPTION( true, std::runtime_error, ss.str().c_str() );
709  }
710  }
711  delete hgradBasis;
712  }
713 }
714 }
715 }
716 
717 #endif
718 
Intrepid2::ArrayTools::scalarMultiplyDataField
static void scalarMultiplyDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties...> outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties...> inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties...> inputFields, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-3, 4, or 5 container inputFields with dimensions (C...
Definition: Intrepid2_ArrayToolsDefScalar.hpp:131
Intrepid2::CellTools::getReferenceEdgeTangent
static void getReferenceEdgeTangent(Kokkos::DynRankView< refEdgeTangentValueType, refEdgeTangentProperties...> refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:391
Intrepid2_DefaultCubatureFactory.hpp
Header file for the abstract base class Intrepid2::DefaultCubatureFactory.
Intrepid2::Impl::OrientationTools::getJacobianOfOrientationMap
static void getJacobianOfOrientationMap(JacobianViewType jacobian, const shards::CellTopology cellTopo, const ordinal_type cellOrt)
Computes jacobian of the parameterization maps of 1- and 2-subcells with orientation.
Intrepid2::OrientationTools::modifyBasisByOrientation
static void modifyBasisByOrientation(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input, const Kokkos::DynRankView< ortValueType, ortProperties...> orts, const BasisPtrType basis)
Modify basis due to orientation.
Definition: Intrepid2_OrientationToolsDefModifyBasis.hpp:133
Intrepid2::FunctionSpaceTools::integrate
static void integrate(Kokkos::DynRankView< outputValueValueType, outputValueProperties...> outputValues, const Kokkos::DynRankView< leftValueValueType, leftValueProperties...> leftValues, const Kokkos::DynRankView< rightValueValueType, rightValueProperties...> rightValues, const bool sumInto=false)
Contracts leftValues and rightValues arrays on the point and possibly space dimensions and stores the...
Definition: Intrepid2_FunctionSpaceToolsDef.hpp:294
Intrepid2::Experimental::ProjectionTools::getHCurlEvaluationPoints
static void getHCurlEvaluationPoints(typename BasisType::scalarViewType evaluationPoints, typename BasisType::scalarViewType curlEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties...> cellOrientations, const BasisType *cellBasis, ProjectionStruct< ExecSpaceType, typename BasisType::scalarType > *projStruct, const EvalPointsType evalPointType=TARGET)
Computes evaluation points for HCurl projection.
Intrepid2::CellTools::mapToReferenceSubcell
static void mapToReferenceSubcell(Kokkos::DynRankView< refSubcellPointValueType, refSubcellPointProperties...> refSubcellPoints, const Kokkos::DynRankView< paramPointValueType, paramPointProperties...> paramPoints, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Computes parameterization maps of 1- and 2-subcells of reference cells.
Definition: Intrepid2_CellToolsDefRefToPhys.hpp:192
Intrepid2_FunctionSpaceTools.hpp
Header file for the Intrepid2::FunctionSpaceTools class.
Intrepid2::RealSpaceTools::clone
static void clone(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input)
Clone input array.
Definition: Intrepid2_RealSpaceToolsDef.hpp:326
Intrepid2::Experimental::ProjectionTools::getHCurlBasisCoeffs
static void getHCurlBasisCoeffs(Kokkos::DynRankView< basisCoeffsValueType, basisCoeffsProperties...> basisCoeffs, const Kokkos::DynRankView< funValsValueType, funValsProperties...> targetAtEvalPoints, const Kokkos::DynRankView< funValsValueType, funValsProperties...> targetCurlAtCurlEvalPoints, const typename BasisType::scalarViewType evaluationPoints, const typename BasisType::scalarViewType curlEvalPoints, const Kokkos::DynRankView< ortValueType, ortProperties...> cellOrientations, const BasisType *cellBasis, ProjectionStruct< ExecSpaceType, typename BasisType::scalarType > *projStruct)
Computes the basis coefficients of the HCurl projection of the target function.
Intrepid2_ArrayTools.hpp
Header file for Intrepid2::ArrayTools class providing utilities for array operations.
Intrepid2::CellTools::getReferenceFaceTangents
static void getReferenceFaceTangents(Kokkos::DynRankView< refFaceTanUValueType, refFaceTanUProperties...> refFaceTanU, Kokkos::DynRankView< refFaceTanVValueType, refFaceTanVProperties...> refFaceTanV, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:428
Intrepid2::Impl::OrientationTools::mapToModifiedReference
static void mapToModifiedReference(outPointViewType outPoints, const refPointViewType refPoints, const shards::CellTopology cellTopo, const ordinal_type cellOrt=0)
Computes modified parameterization maps of 1- and 2-subcells with orientation.
Definition: Intrepid2_OrientationToolsDefModifyPoints.hpp:254
Intrepid2::CellTools::getReferenceFaceNormal
static void getReferenceFaceNormal(Kokkos::DynRankView< refFaceNormalValueType, refFaceNormalProperties...> refFaceNormal, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to faces of 3D reference cell.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:504
Generated by   doxygen 1.8.5