From b6a1cb13afe3405ddc3ab53bafa80357c4b5c0fe Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Robert=20Kl=C3=B6fkorn?= <robertk@dune-project.org>
Date: Mon, 3 Apr 2006 17:48:17 +0000
Subject: [PATCH] also removed operator, because feoperator.hh not existing
anymore.
[[Imported from SVN: r4439]]
---
fem/Makefile.am | 2 +-
fem/operator/.gitignore | 3 -
fem/operator/Makefile.am | 6 -
fem/operator/laplace.hh | 204 -----------------------
fem/operator/massmatrix.hh | 323 -------------------------------------
5 files changed, 1 insertion(+), 537 deletions(-)
delete mode 100644 fem/operator/.gitignore
delete mode 100644 fem/operator/Makefile.am
delete mode 100644 fem/operator/laplace.hh
delete mode 100644 fem/operator/massmatrix.hh
diff --git a/fem/Makefile.am b/fem/Makefile.am
index 6c18cd541..01252021e 100644
--- a/fem/Makefile.am
+++ b/fem/Makefile.am
@@ -1,6 +1,6 @@
# $Id$
-SUBDIRS = norms operator transfer
+SUBDIRS = norms transfer
femincludedir = $(includedir)/dune/fem
feminclude_HEADERS = l2projection.hh
diff --git a/fem/operator/.gitignore b/fem/operator/.gitignore
deleted file mode 100644
index 55880268f..000000000
--- a/fem/operator/.gitignore
+++ /dev/null
@@ -1,3 +0,0 @@
-Makefile
-Makefile.in
-semantic.cache
\ No newline at end of file
diff --git a/fem/operator/Makefile.am b/fem/operator/Makefile.am
deleted file mode 100644
index c9380eb89..000000000
--- a/fem/operator/Makefile.am
+++ /dev/null
@@ -1,6 +0,0 @@
-# $Id$
-
-femopdir = $(includedir)/dune/fem/operator
-femop_HEADERS = laplace.hh massmatrix.hh
-
-include $(top_srcdir)/am/global-rules
diff --git a/fem/operator/laplace.hh b/fem/operator/laplace.hh
deleted file mode 100644
index a7c499ca0..000000000
--- a/fem/operator/laplace.hh
+++ /dev/null
@@ -1,204 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_LAPLACE_HH
-#define DUNE_LAPLACE_HH
-
-#include <dune/fem/feoperator.hh>
-#include <dune/fem/feop/spmatrix.hh>
-
-#include <dune/quadrature/quadraturerules.hh>
-
-namespace Dune
-{
-
- /** \brief The Laplace operator
- */
- template <class DiscFunctionType, class TensorType, int polOrd>
- class LaplaceFEOp :
- public FiniteElementOperator<DiscFunctionType,
- SparseRowMatrix<double>,
- LaplaceFEOp<DiscFunctionType,TensorType, polOrd> > {
-
- //! The corresponding function space type
- typedef typename DiscFunctionType::FunctionSpaceType FunctionSpaceType;
-
- //! The grid
- typedef typename FunctionSpaceType::GridType GridType;
-
- //! The grid's dimension
- enum { dim = GridType::dimension };
-
- //! ???
- typedef typename FunctionSpaceType::JacobianRangeType JacobianRange;
-
- //! ???
- typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
-
- //! ???
- typedef typename FiniteElementOperator<DiscFunctionType,
- SparseRowMatrix<double>,
- LaplaceFEOp<DiscFunctionType, TensorType, polOrd> >::OpMode OpMode;
-
-
- mutable JacobianRange grad;
- mutable JacobianRange othGrad;
- mutable JacobianRange mygrad[4];
-
- public:
-
- //! ???
- DiscFunctionType *stiffFunktion_;
-
- //! ???
- TensorType *stiffTensor_;
-
- //! ???
- LaplaceFEOp( const typename DiscFunctionType::FunctionSpaceType &f, OpMode opMode ) :
- FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,LaplaceFEOp<DiscFunctionType,TensorType, polOrd> >( f, opMode ) ,
- stiffFunktion_(NULL), stiffTensor_(NULL)
- {}
-
- //! ???
- LaplaceFEOp( const DiscFunctionType &stiff, const typename DiscFunctionType::FunctionSpaceType &f, OpMode opMode ) :
- FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,LaplaceFEOp<DiscFunctionType,TensorType, polOrd> >( f, opMode ) ,
- stiffFunktion_(&stiff), stiffTensor_(NULL)
- {}
-
- //! ???
- LaplaceFEOp( TensorType &stiff, const typename DiscFunctionType::FunctionSpaceType &f, OpMode opMode ) :
- FiniteElementOperator<DiscFunctionType,SparseRowMatrix<double>,LaplaceFEOp<DiscFunctionType,TensorType, polOrd> >( f, opMode ) ,
- stiffFunktion_(NULL), stiffTensor_(&stiff)
- {}
-
- //! Returns the actual matrix if it is assembled
- const SparseRowMatrix<double>* getMatrix() const {
- assert(this->matrix_);
- return this->matrix_;
- }
-
- //! Creates a new empty matrix
- SparseRowMatrix<double>* newEmptyMatrix( ) const
- {
- return new SparseRowMatrix<double>( this->functionSpace_.size () ,
- this->functionSpace_.size () ,
- 15 * dim);
- }
-
- //! Prepares the local operator before calling apply()
- void prepareGlobal ( const DiscFunctionType &arg, DiscFunctionType &dest )
- {
- this->arg_ = &arg.argument();
- this->dest_ = &dest.destination();
- assert(this->arg_ != 0); assert(this->dest_ != 0);
- dest.clear();
- }
-
- //! ???
- template <class EntityType>
- double getLocalMatrixEntry( EntityType &entity, const int i, const int j ) const
- {
- enum { dim = GridType::dimension };
- typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
-
- const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
- // Get quadrature rule
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(entity.geometry().type(), polOrd);
-
- double val = 0.;
- for ( int pt=0; pt < quad.size(); pt++ )
- {
- baseSet.jacobian(i,quad[pt].position(),grad);
-
- // calc Jacobian inverse before volume is evaluated
- const FieldMatrix<double,dim,dim>& inv = entity.geometry().jacobianInverse(quad[pt].position());
- const double vol = entity.geometry().integrationElement(quad[pt].position());
-
- // multiply with transpose of jacobian inverse
- JacobianRange tmp(0);
- inv.umv(grad[0], tmp[0]);
- grad[0] = tmp[0];
-
- if( i != j )
- {
- baseSet.jacobian(j,quad[pt].position(),othGrad);
-
- // multiply with transpose of jacobian inverse
- JacobianRange tmp(0);
- inv.umv(othGrad[0], tmp[0]);
- othGrad[0] = tmp[0];
-
- val += ( grad[0] * othGrad[0] ) * quad[pt].weight() * vol;
- }
- else
- {
- val += ( grad[0] * grad[0] ) * quad[pt].weight() * vol;
- }
- }
-
- return val;
- }
-
- //! ???
- template < class EntityType, class MatrixType>
- void getLocalMatrix( EntityType &entity, const int matSize, MatrixType& mat) const {
- enum { dim = GridType::dimension };
-
- typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
-
- const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
-
- int i,j;
-
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- mat[j][i]=0.0;
-
- // Get quadrature rule
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(entity.geometry().type(), polOrd);
-
- for ( int pt=0; pt < quad.size(); pt++ ) {
-
- // calc Jacobian inverse before volume is evaluated
- const FieldMatrix<double,dim,dim>& inv = entity.geometry().jacobianInverse(quad[pt].position());
- const double vol = entity.geometry().integrationElement(quad[pt].position());
-
- for(i=0; i<matSize; i++) {
-
- baseSet.jacobian(i,quad[pt].position(),mygrad[i]);
-
- // multiply with transpose of jacobian inverse
- JacobianRange tmp(0);
- inv.umv(mygrad[i][0], tmp[0]);
- mygrad[i][0] = tmp[0];
- }
-
- typename FunctionSpaceType::RangeType ret;
-
- if(stiffTensor_) {
- stiffTensor_->evaluate(entity.geometry().global(quad[pt].position()),ret);
- ret[0] *= quad[pt].weight();
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- mat[j][i] += ( mygrad[i][0] * mygrad[j][0] ) * ret[0] * vol;
- }
- else{
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- mat[j][i] += ( mygrad[i][0] * mygrad[j][0] ) * quad[pt].weight() * vol;
- }
-
-
-
- }
-
- for(i=0; i<matSize; i++)
- for (j=matSize-1; j>i; j--)
- mat[j][i] = mat[i][j];
-
- }
-
- }; // end class
-
-} // end namespace
-
-#endif
diff --git a/fem/operator/massmatrix.hh b/fem/operator/massmatrix.hh
deleted file mode 100644
index a260ba94c..000000000
--- a/fem/operator/massmatrix.hh
+++ /dev/null
@@ -1,323 +0,0 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-#ifndef DUNE_MASSMATRIX_HH
-#define DUNE_MASSMATRIX_HH
-
-#include <dune/fem/feoperator.hh>
-#include <dune/fem/feop/spmatrix.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/matrixindexset.hh>
-
-#include <dune/quadrature/quadraturerules.hh>
-
-namespace Dune {
-
- /** \brief The mass matrix
- *
- * \tparam polOrd The quadrature order
- *
- * \todo Giving the quadrature order as a template parameter is
- * a hack. It would be better to determine the optimal order automatically.
- */
- template <class DiscFunctionType, int polOrd>
- class MassMatrixFEOp :
-
- public FiniteElementOperator<DiscFunctionType,
- SparseRowMatrix<double>,
- MassMatrixFEOp<DiscFunctionType, polOrd> > {
- typedef typename DiscFunctionType::FunctionSpaceType FunctionSpaceType;
- typedef typename FiniteElementOperator<DiscFunctionType,
- SparseRowMatrix<double>,
- MassMatrixFEOp<DiscFunctionType, polOrd> >::OpMode OpMode;
-
- typedef typename FunctionSpaceType::RangeField RangeFieldType;
- typedef typename FunctionSpaceType::GridType GridType;
- typedef typename FunctionSpaceType::Range RangeType;
-
- public:
-
- //! Returns the actual matrix if it is assembled
- /** \todo Should this be in a base class? */
- const SparseRowMatrix<double>* getMatrix() const {
- assert(this->matrix_);
- return this->matrix_;
- }
-
- //! Constructor
- MassMatrixFEOp( const typename DiscFunctionType::FunctionSpace &f, OpMode opMode ) : //= ON_THE_FLY ) :
- FiniteElementOperator<DiscFunctionType,
- SparseRowMatrix<double>,
- MassMatrixFEOp<DiscFunctionType, polOrd> >( f, opMode )
- {}
-
- //! ???
- SparseRowMatrix<double>* newEmptyMatrix( ) const {
- return new SparseRowMatrix<double>( this->functionSpace_.size () ,
- this->functionSpace_.size () ,
- 30);
- }
-
- //! ???
- template <class EntityType>
- double getLocalMatrixEntry( EntityType &entity, const int i, const int j ) const
- {
- typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
- const int dim = GridType::dimension;
- double val = 0;
-
- FieldVector<double, dim> tmp(1.0);
-
- const BaseFunctionSetType & baseSet
- = this->functionSpace_.getBaseFunctionSet( entity );
-
- RangeType v1 (0.0);
- RangeType v2 (0.0);
-
- // Get quadrature rule
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(entity.geometry().type(), polOrd);
-
- const double vol = entity.geometry().integrationElement( tmp );
- for ( int pt=0; pt < quad.size(); pt++ )
- {
- baseSet.eval(i,quad,pt,v1);
- baseSet.eval(j,quad,pt,v2);
- val += ( v1 * v2 ) * quad[pt].weight();
- }
-
- return val*vol;
- }
-
- //! ???
- template < class EntityType, class MatrixType>
- void getLocalMatrix( EntityType &entity, const int matSize, MatrixType& mat) const
- {
- int i,j;
- typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
- const int dim = GridType::dimension;
- const BaseFunctionSetType & baseSet
- = this->functionSpace_.getBaseFunctionSet( entity );
-
- /** \todo What's the correct type here? */
- static FieldVector<double, dim> tmp(1.0);
- const double vol = entity.geometry().integrationElement(tmp);
-
- static RangeType v[30];
- // Check magic constant. Otherwise program will fail in loop below
- assert(matSize <= 30);
-
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- mat[j][i]=0.0;
-
- // Get quadrature rule
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(entity.geometry().type(), polOrd);
-
- for ( int pt=0; pt < quad.size(); pt++ )
- {
- for(i=0; i<matSize; i++)
- baseSet.eval(i,quad,pt,v[i]);
-
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- mat[j][i] += ( v[i] * v[j] ) * quad[pt].weight();
- }
-
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- mat[j][i] *= vol;
-
- for(i=0; i<matSize; i++)
- for (j=matSize; j>i; j--)
- mat[j][i] = mat[i][j];
-
- // * for debugging only
- // for (int i = 0; i < matSize; ++i) {
- //assert(mat[i][i] > 0.0039 && mat[i][i] < 0.004);
- //}
- }
-
- protected:
- DiscFunctionType *_tmp;
- };
-
-
-
-
-
-
- /** \brief The mass matrix operator for block matrices
- *
- * \tparam polOrd The quadrature order
- *
- * This is an implementation of the mass matrix operator using the ISTL-classes
- *
- * \todo Giving the quadrature order as a template parameter is
- * a hack. It would be better to determine the optimal order automatically.
- */
- template <class GridType, int blocksize, int polOrd>
- class MassMatrix {
-
- enum {dim = GridType::dimension};
-
- enum {elementOrder = 1};
-
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
- //!
- typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
- public:
-
- /** \todo Does actually belong into the base class */
- BCRSMatrix<MatrixBlock>* matrix_;
-
- /** \todo Does actually belong into the base class */
- const GridType* grid_;
-
- //! ???
- MassMatrix(const GridType &grid) :
- grid_(&grid)
- {}
-
- //! Returns the actual matrix if it is assembled
- const BCRSMatrix<MatrixBlock>* getMatrix() const {
- assert(this->matrix_);
- return this->matrix_;
- }
-
- void getNeighborsPerVertex(MatrixIndexSet& nb) {
-
- int i, j;
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
- int n = indexSet.size(dim);
-
- nb.resize(n, n);
-
- typedef typename GridType::template Codim<0>::LevelIterator LevelIterator;
- LevelIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
- LevelIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
-
- for (; it!=endit; ++it) {
-
- for (i=0; i<it->template count<dim>(); i++) {
-
- for (j=0; j<it->template count<dim>(); j++) {
-
- int iIdx = indexSet.template subIndex<dim>(*it, i);
- int jIdx = indexSet.template subIndex<dim>(*it, j);
-
- nb.add(iIdx, jIdx);
-
- }
-
- }
-
- }
-
- }
-
- void assembleMatrix() {
-
- typedef typename GridType::template Codim<0>::LevelIterator LevelIterator;
-
- const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
- int n = indexSet.size(GridType::dimension);
-
- MatrixIndexSet neighborsPerVertex;
- getNeighborsPerVertex(neighborsPerVertex);
-
- matrix_ = new BCRSMatrix<MatrixBlock>(n, n, BCRSMatrix<MatrixBlock>::random);
-
- neighborsPerVertex.exportIdx(*matrix_);
- (*matrix_) = 0;
-
- LevelIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
- LevelIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
- enum {maxnumOfBaseFct = 30};
-
- Matrix<MatrixBlock> mat;
-
- for( ; it != endit; ++it ) {
-
- // const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
- // const int numOfBaseFct = baseSet.numBaseFunctions();
- const LagrangeShapeFunctionSet<double, double, dim> & baseSet
- = Dune::LagrangeShapeFunctions<double, double, dim>::general(it->geometry().type(), elementOrder);
-
- const int numOfBaseFct = baseSet.size();
- //printf("%d base functions on element\n", numOfBaseFct);
- mat.resize(numOfBaseFct, numOfBaseFct);
- // setup matrix
- getLocalMatrix( *it, numOfBaseFct, mat);
-
- for(int i=0; i<numOfBaseFct; i++) {
-
- //int row = functionSpace_.mapToGlobal( *it , i );
- int row = indexSet.template subIndex<dim>(*it,i);
- for (int j=0; j<numOfBaseFct; j++ ) {
-
- //int col = functionSpace_.mapToGlobal( *it , j );
- int col = indexSet.template subIndex<dim>(*it,j);
- (*matrix_)[row][col] += mat[i][j];
-
- }
- }
- }
-
- }
-
- //! ???
- template < class EntityType, class MatrixType>
- void getLocalMatrix( EntityType &entity, const int matSize, MatrixType& mat) const
- {
- int i,j;
- const LagrangeShapeFunctionSet<double, double, dim> & baseSet
- = Dune::LagrangeShapeFunctions<double, double, dim>::general(entity.geometry().type(), elementOrder);
- // Clear local scalar matrix
- Matrix<double> scalarMat(matSize, matSize);
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++ )
- scalarMat[i][j]=0.0;
-
- // Get quadrature rule
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(entity.geometry().type(), polOrd);
-
- for (unsigned int pt=0; pt < quad.size(); pt++ ) {
-
- // Get position of the quadrature point
- const FieldVector<double,dim>& quadPos = quad[pt].position();
-
- // The factor in the integral transformation formula
- const double integrationElement = entity.geometry().integrationElement(quadPos);
-
- double v[matSize];
- for(i=0; i<matSize; i++)
- v[i] = baseSet[i].evaluateFunction(0,quadPos);
-
- for(i=0; i<matSize; i++)
- for (int j=0; j<=i; j++ )
- scalarMat[i][j] += ( v[i] * v[j] ) * quad[pt].weight() * integrationElement;
- }
-
- // Clear local matrix
- for(i=0; i<matSize; i++)
- for (j=0; j<matSize; j++)
- mat[i][j] = 0;
-
- // Turn scalar triangular matrix into symmetric block matrix
- for(i=0; i<matSize; i++)
- for (j=0; j<=i; j++)
- for (int k=0; k<blocksize; k++)
- mat[i][j][k][k] = mat[j][i][k][k] = scalarMat[i][j];
-
- }
-
-
- }; // end class
-
-
-
-
-} // namespace Dune
-
-#endif
--
GitLab