Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
00030 
00031 #ifndef SUNDANCE_MAXIMALQUADRATUREINTEGRAL_H
00032 #define SUNDANCE_MAXIMALQUADRATUREINTEGRAL_H
00033 
00034 #include "SundanceDefs.hpp"
00035 #include "SundanceQuadratureIntegralBase.hpp"
00036 
00037 namespace Sundance
00038 {
00039 
00040 using namespace Teuchos;
00041 
00042 /** 
00043  * Integration by quadrature over a maximal cell  
00044  */
00045 class MaximalQuadratureIntegral : public ElementIntegral
00046 {
00047 public:
00048   /** Construct a zero-form to be computed by quadrature */
00049   MaximalQuadratureIntegral(
00050     const CellType& maxCellType,
00051     const QuadratureFamily& quad,
00052     const ParametrizedCurve& globalCurve,
00053     const Mesh& mesh,
00054     int verb);
00055 
00056   /** Construct a one form to be computed by quadrature */
00057   MaximalQuadratureIntegral(
00058     const CellType& maxCellType,
00059     const BasisFamily& testBasis,
00060     int alpha,
00061     int testDerivOrder,
00062     const QuadratureFamily& quad,
00063     const ParametrizedCurve& globalCurve,
00064     const Mesh& mesh,
00065     int verb);
00066 
00067   /** Construct a two-form to be computed by quadrature */
00068   MaximalQuadratureIntegral(
00069     const CellType& maxCellType,
00070     const BasisFamily& testBasis,
00071     int alpha,
00072     int testDerivOrder,
00073     const BasisFamily& unkBasis,
00074     int beta,
00075     int unkDerivOrder,
00076     const QuadratureFamily& quad,
00077     const ParametrizedCurve& globalCurve,
00078     const Mesh& mesh,
00079     int verb);
00080 
00081   /** virtual dtor */
00082   virtual ~MaximalQuadratureIntegral(){;}
00083       
00084      /** */
00085   virtual void transform(const CellJacobianBatch& JTrans,
00086     const CellJacobianBatch& JVol,
00087     const Array<int>& isLocalFlag,
00088     const Array<int>& facetNum,
00089     const RCP<Array<int> >& cellLIDs,
00090     const double* const coeff,
00091     RCP<Array<double> >& A) const 
00092     {
00093       if (order()==2) transformTwoForm(JTrans, JVol, facetNum, cellLIDs,coeff, A);
00094       else if (order()==1) transformOneForm(JTrans, JVol, facetNum, cellLIDs,coeff, A);
00095       else transformZeroForm(JTrans, JVol, isLocalFlag, facetNum, cellLIDs,coeff, A);
00096     }
00097 
00098   /** */
00099   virtual void transformZeroForm(const CellJacobianBatch& JTrans,
00100     const CellJacobianBatch& JVol,
00101     const Array<int>& isLocalFlag,
00102     const Array<int>& facetIndex,
00103     const RCP<Array<int> >& cellLIDs,
00104     const double* const coeff,
00105     RCP<Array<double> >& A) const ;
00106   
00107   /** */
00108   virtual void transformTwoForm(const CellJacobianBatch& JTrans,
00109     const CellJacobianBatch& JVol,
00110     const Array<int>& facetIndex,
00111     const RCP<Array<int> >& cellLIDs,
00112     const double* const coeff,
00113     RCP<Array<double> >& A) const ;
00114   
00115   /** */
00116   void transformOneForm(const CellJacobianBatch& JTrans,
00117     const CellJacobianBatch& JVol,
00118     const Array<int>& facetIndex,
00119     const RCP<Array<int> >& cellLIDs,
00120     const double* const coeff,
00121     RCP<Array<double> >& A) const ;
00122   
00123 private:
00124 
00125   /** Do the integration by summing reference quantities over quadrature
00126    * points and then transforming the sum to physical quantities.  */
00127   void transformSummingFirst(int nCells,
00128     const Array<int>& facetIndex,
00129     const RCP<Array<int> >& cellLIDs,
00130     const double* const GPtr,
00131     const double* const coeff,
00132     RCP<Array<double> >& A) const ;
00133 
00134   /** Do the integration by transforming to physical coordinates 
00135    * at each quadrature point, and then summing */
00136   void transformSummingLast(int nCells,
00137     const Array<int>& facetIndex,
00138     const RCP<Array<int> >& cellLIDs,
00139     const double* const GPtr,
00140     const double* const coeff,
00141     RCP<Array<double> >& A) const ;
00142 
00143   /** Determine whether to do this batch of integrals using the
00144    * sum-first method or the sum-last method */
00145   bool useSumFirstMethod() const {return useSumFirstMethod_;}
00146       
00147   /** */
00148   inline double& wValue(int q, int testDerivDir, int testNode,
00149     int unkDerivDir, int unkNode)
00150     {return W_[unkNode
00151         + nNodesUnk()
00152         *(testNode + nNodesTest()
00153           *(unkDerivDir + nRefDerivUnk()
00154             *(testDerivDir + nRefDerivTest()*q)))];}
00155 
00156       
00157 
00158   /** */
00159   inline const double& wValue(int facetCase, 
00160     int q, 
00161     int testDerivDir, int testNode,
00162     int unkDerivDir, int unkNode) const 
00163     {
00164       return W_[unkNode
00165         + nNodesUnk()
00166         *(testNode + nNodesTest()
00167           *(unkDerivDir + nRefDerivUnk()
00168             *(testDerivDir + nRefDerivTest()*q)))];
00169     }
00170       
00171   /** */
00172   inline double& wValue(int q, int testDerivDir, int testNode)
00173     {return W_[testNode + nNodesTest()*(testDerivDir + nRefDerivTest()*q)];}
00174 
00175 
00176   /** */
00177   inline const double& wValue(int q, int testDerivDir, int testNode) const 
00178     {return W_[testNode + nNodesTest()*(testDerivDir + nRefDerivTest()*q)];}
00179 
00180   /* */
00181   Array<double> W_;
00182 
00183   /* */
00184   bool useSumFirstMethod_;
00185       
00186   /** For ACI (ACI = Adaptive Cell Integration), store the reference integral values for one form
00187    * The indexes facet, quadPoints, nRefDerivTest , nNodesTest */
00188   Array<Array<Array<double> > >  W_ACI_F1_;
00189 
00190   /** For ACI (ACI = Adaptive Cell Integration), store the reference integral values for two form
00191    * The indexes quadPoints, nRefDerivTest , nNodesTest , nRefDerivUnk , nNodesUnk */
00192   Array<Array<Array<Array<Array<double> > > > > W_ACI_F2_;
00193 
00194   /** The quadrature family needed for special integration (ACI)*/
00195   QuadratureFamily quad_;
00196 
00197   /** The quadrature points*/
00198   Array<Point> quadPts_;
00199 
00200   /** The standard weights (in case of ACI these might change)*/
00201   Array<double> quadWeights_;
00202 };
00203 }
00204 
00205 
00206 #endif
SundanceMaximalQuadratureIntegral.hpp
