From 2fd3339195d1c1f77dce30a211b612defb2a50d9 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 30 Nov 2005 18:03:25 +0000
Subject: [PATCH] using new base function implementations
[[Imported from SVN: r597]]
---
src/rodassembler.cc | 199 +++++++++++++++++++++++++-------------------
src/rodassembler.hh | 25 ++----
staticrod.cc | 52 ++++--------
3 files changed, 133 insertions(+), 143 deletions(-)
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 020cc3d6..c8ec0113 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -5,21 +5,22 @@
#include <dune/quadrature/quadraturerules.hh>
+#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
-
-template <class FunctionSpaceType, int polOrd>
-void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+template <class GridType, int polOrd>
+void Dune::RodAssembler<GridType, polOrd>::
getNeighborsPerVertex(MatrixIndexSet& nb) const
{
const int gridDim = GridType::dimension;
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
int i, j;
- int n = grid_->size(grid_->maxlevel(), gridDim);
+ int n = grid_->size(grid_->maxLevel(), gridDim);
nb.resize(n, n);
- ElementIterator it = grid_->template lbegin<0>( grid_->maxlevel() );
- ElementIterator endit = grid_->template lend<0> ( grid_->maxlevel() );
+ ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
+ ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
for (; it!=endit; ++it) {
@@ -27,8 +28,10 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const
for (j=0; j<it->template count<gridDim>(); j++) {
- int iIdx = it->template subIndex<gridDim>(i);
- int jIdx = it->template subIndex<gridDim>(j);
+ //int iIdx = it->template subIndex<gridDim>(i);
+ //int jIdx = it->template subIndex<gridDim>(j);
+ int iIdx = indexSet.template subIndex<gridDim>(*it,i);
+ int jIdx = indexSet.template subIndex<gridDim>(*it,j);
nb.add(iIdx, jIdx);
@@ -40,26 +43,29 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const
}
-template <class FunctionSpaceType, int polOrd>
-void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+template <class GridType, int polOrd>
+void Dune::RodAssembler<GridType, polOrd>::
assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
BCRSMatrix<MatrixBlock>& matrix)
{
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
matrix = 0;
- ElementIterator it = grid_->template lbegin<0>( grid_->maxlevel() );
- ElementIterator endit = grid_->template lend<0> ( grid_->maxlevel() );
-
+ ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
+ ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
+
Matrix<MatrixBlock> mat;
for( ; it != endit; ++it ) {
- const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
- const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();
+ //const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
+ const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->geometry().type(), elementOrder);
+ const int numOfBaseFct = baseSet.size();
mat.resize(numOfBaseFct, numOfBaseFct);
@@ -67,7 +73,8 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
BlockVector<FieldVector<double, blocksize> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+ //localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+ localSolution[i] = sol[indexSet.template subIndex<gridDim>(*it,i)];
// setup matrix
getLocalMatrix( *it, localSolution, numOfBaseFct, mat);
@@ -75,11 +82,13 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
- int row = functionSpace_.mapToGlobal( *it , i );
+ //int row = functionSpace_.mapToGlobal( *it , i );
+ int row = indexSet.template subIndex<gridDim>(*it,i);
for (int j=0; j<numOfBaseFct; j++ ) {
- int col = functionSpace_.mapToGlobal( *it , j );
+ //int col = functionSpace_.mapToGlobal( *it , j );
+ int col = indexSet.template subIndex<gridDim>(*it,j);
matrix[row][col] += mat[i][j];
}
@@ -94,14 +103,15 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
-template <class FunctionSpaceType, int polOrd>
+template <class GridType, int polOrd>
template <class MatrixType>
-void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+void Dune::RodAssembler<GridType, polOrd>::
getLocalMatrix( EntityType &entity,
const BlockVector<FieldVector<double, blocksize> >& localSolution,
const int matSize, MatrixType& localMat) const
{
-
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
+
/* ndof is the number of vectors of the element */
int ndof = matSize;
@@ -109,9 +119,10 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<matSize; j++)
localMat[i][j] = 0;
- typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
- const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
-
+ //const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
+ const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(entity.geometry().type(), elementOrder);
+
// Get quadrature rule
const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(entity.geometry().type(), polOrd);
@@ -122,7 +133,7 @@ getLocalMatrix( EntityType &entity,
const FieldVector<double,gridDim>& quadPos = quad[ip].position();
// calc Jacobian inverse before integration element is evaluated
- const FieldMatrix<double,gridDim,gridDim>& inv = entity.geometry().jacobianInverse(quadPos);
+ const FieldMatrix<double,gridDim,gridDim>& inv = entity.geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = entity.geometry().integrationElement(quadPos);
/* Compute the weight of the current integration point */
@@ -131,31 +142,33 @@ getLocalMatrix( EntityType &entity,
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
- Array<JacobianRange> shapeGrad(ndof);
+ Array<FieldVector<double,gridDim> > shapeGrad(ndof);
for (int dof=0; dof<ndof; dof++) {
- baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+ //baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+ for (int i=0; i<gridDim; i++)
+ shapeGrad[dof][i] = baseSet[dof].evaluateDerivative(0,i,quadPos);
// multiply with jacobian inverse
- JacobianRange tmp(0);
- inv.umv(shapeGrad[dof][0], tmp[0]);
- shapeGrad[dof][0] = tmp[0];
+ FieldVector<double,gridDim> tmp(0);
+ inv.umv(shapeGrad[dof], tmp);
+ shapeGrad[dof] = tmp;
}
- RangeType shapeFunction[matSize];
+ double shapeFunction[matSize];
for(int i=0; i<matSize; i++)
- baseSet.eval(i,quadPos,shapeFunction[i]);
-
+ //baseSet.eval(i,quadPos,shapeFunction[i]);
+ shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
+
// //////////////////////////////////
// Interpolate
// //////////////////////////////////
- double x_s = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*shapeGrad[1][0][0];
- double y_s = localSolution[0][1]*shapeGrad[0][0][0] + localSolution[1][1]*shapeGrad[1][0][0];
- //double theta_s = localSolution[0][2]*shapeGrad[0][0][0] + localSolution[1][2]*shapeGrad[1][0][0];
+ double x_s = localSolution[0][0]*shapeGrad[0][0] + localSolution[1][0]*shapeGrad[1][0];
+ double y_s = localSolution[0][1]*shapeGrad[0][0] + localSolution[1][1]*shapeGrad[1][0];
double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1];
@@ -164,51 +177,51 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<matSize; j++) {
// \partial J^2 / \partial x_i \partial x_j
- localMat[i][j][0][0] += weight * shapeGrad[i][0][0] * shapeGrad[j][0][0]
+ localMat[i][j][0][0] += weight * shapeGrad[i][0] * shapeGrad[j][0]
* (A1 * cos(theta) * cos(theta) + A3 * sin(theta) * sin(theta));
// \partial J^2 / \partial x_i \partial y_j
- localMat[i][j][0][1] += weight * shapeGrad[i][0][0] * shapeGrad[j][0][0]
+ localMat[i][j][0][1] += weight * shapeGrad[i][0] * shapeGrad[j][0]
* (-A1 + A3) * sin(theta)* cos(theta);
// \partial J^2 / \partial x_i \partial theta_j
- localMat[i][j][0][2] += weight * shapeGrad[i] * shapeFunction[j]
+ localMat[i][j][0][2] += weight * shapeGrad[i][0] * shapeFunction[j]
* (-A1 * (x_s*sin(theta) + y_s*cos(theta)) * cos(theta)
- A1* (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
+A3 * (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
+A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta));
// \partial J^2 / \partial y_i \partial x_j
- localMat[i][j][1][0] += weight * shapeGrad[i][0][0] * shapeGrad[j][0][0]
+ localMat[i][j][1][0] += weight * shapeGrad[i][0] * shapeGrad[j][0]
* (-A1 * sin(theta)* cos(theta) + A3 * cos(theta) * sin(theta));
// \partial J^2 / \partial y_i \partial y_j
- localMat[i][j][1][1] += weight * shapeGrad[i] * shapeGrad[j]
+ localMat[i][j][1][1] += weight * shapeGrad[i][0] * shapeGrad[j][0]
* (A1 * sin(theta)*sin(theta) + A3 * cos(theta)*cos(theta));
// \partial J^2 / \partial y_i \partial theta_j
- localMat[i][j][1][2] += weight * shapeGrad[i] * shapeFunction[j]
+ localMat[i][j][1][2] += weight * shapeGrad[i][0] * shapeFunction[j]
* (A1 * (x_s * sin(theta) + y_s * cos(theta)) * sin(theta)
-A1 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
+A3 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
-A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta));
// \partial J^2 / \partial theta_i \partial x_j
- localMat[i][j][2][0] += weight * shapeFunction[i] * shapeGrad[j]
+ localMat[i][j][2][0] += weight * shapeFunction[i] * shapeGrad[j][0]
* (-A1 * (x_s*sin(theta) + y_s*cos(theta)) * cos(theta)
- A1* (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
+A3 * (x_s*cos(theta) - y_s*sin(theta)) * sin(theta)
+A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta));
// \partial J^2 / \partial theta_i \partial y_j
- localMat[i][j][2][1] += weight * shapeFunction[i] * shapeGrad[j]
+ localMat[i][j][2][1] += weight * shapeFunction[i] * shapeGrad[j][0]
* (A1 * (x_s * sin(theta) + y_s * cos(theta)) * sin(theta)
-A1 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
+A3 * (x_s * cos(theta) - y_s * sin(theta)) * cos(theta)
-A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta));
// \partial J^2 / \partial theta_i \partial theta_j
- localMat[i][j][2][2] += weight * B * shapeGrad[i] * shapeGrad[j];
+ localMat[i][j][2][2] += weight * B * shapeGrad[i][0] * shapeGrad[j][0];
localMat[i][j][2][2] += weight * shapeFunction[i] * shapeFunction[j]
* (+ A1 * (x_s*sin(theta) + y_s*cos(theta)) * (x_s*sin(theta) + y_s*cos(theta))
+ A1 * (x_s*cos(theta) - y_s*sin(theta)) * (-x_s*cos(theta)+ y_s*sin(theta))
@@ -252,12 +265,13 @@ getLocalMatrix( EntityType &entity,
}
-template <class FunctionSpaceType, int polOrd>
-void Dune::RodAssembler<FunctionSpaceType, polOrd>::
+template <class GridType, int polOrd>
+void Dune::RodAssembler<GridType, polOrd>::
assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
BlockVector<FieldVector<double, blocksize> >& grad) const
{
- const int maxlevel = grid_->maxlevel();
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
+ const int maxlevel = grid_->maxLevel();
if (sol.size()!=grid_->size(maxlevel, gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
@@ -272,13 +286,16 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
for (; it!=endIt; ++it) {
// Extract local solution on this element
- const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
- const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();
+ //const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
+ const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->geometry().type(), elementOrder);
+ const int numOfBaseFct = baseSet.size();
FieldVector<double, blocksize> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+ //localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+ localSolution[i] = sol[indexSet.template subIndex<gridDim>(*it,i)];
// Get quadrature rule
const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().type(), polOrd);
@@ -288,7 +305,7 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
// Local position of the quadrature point
const FieldVector<double,gridDim>& quadPos = quad[pt].position();
- const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverse(quadPos);
+ const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = it->geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement;
@@ -296,31 +313,33 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
// ///////////////////////////////////////
// Compute deformation gradient
// ///////////////////////////////////////
- Array<JacobianRange> shapeGrad(numOfBaseFct);
+ Array<FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
for (int dof=0; dof<numOfBaseFct; dof++) {
- baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
-
+ //baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+ for (int i=0; i<gridDim; i++)
+ shapeGrad[dof][i] = baseSet[dof].evaluateDerivative(0,i,quadPos);
+
// multiply with jacobian inverse
- JacobianRange tmp(0);
- inv.umv(shapeGrad[dof][0], tmp[0]);
- shapeGrad[dof][0] = tmp[0];
+ FieldVector<double,gridDim> tmp(0);
+ inv.umv(shapeGrad[dof], tmp);
+ shapeGrad[dof] = tmp;
}
// Get the value of the shape functions
- RangeType shapeFunction[2];
+ double shapeFunction[2];
for(int i=0; i<2; i++)
- baseSet.eval(i,quadPos,shapeFunction[i]);
+ shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
// //////////////////////////////////
// Interpolate
// //////////////////////////////////
- double x_s = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*shapeGrad[1][0][0];
- double y_s = localSolution[0][1]*shapeGrad[0][0][0] + localSolution[1][1]*shapeGrad[1][0][0];
- double theta_s = localSolution[0][2]*shapeGrad[0][0][0] + localSolution[1][2]*shapeGrad[1][0][0];
+ double x_s = localSolution[0][0]*shapeGrad[0][0] + localSolution[1][0]*shapeGrad[1][0];
+ double y_s = localSolution[0][1]*shapeGrad[0][0] + localSolution[1][1]*shapeGrad[1][0];
+ double theta_s = localSolution[0][2]*shapeGrad[0][0] + localSolution[1][2]*shapeGrad[1][0];
double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1];
@@ -333,18 +352,19 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
for (int dof=0; dof<numOfBaseFct; dof++) {
- int globalDof = functionSpace_.mapToGlobal(*it,dof);
+ //int globalDof = functionSpace_.mapToGlobal(*it,dof);
+ int globalDof = indexSet.template subIndex<gridDim>(*it,dof);
//printf("globalDof: %d partA1: %g partA3: %g\n", globalDof, partA1, partA3);
// \partial J / \partial x^i
- grad[globalDof][0] += weight * (partA1 * cos(theta) + partA3 * sin(theta)) * shapeGrad[dof][0][0];
+ grad[globalDof][0] += weight * (partA1 * cos(theta) + partA3 * sin(theta)) * shapeGrad[dof][0];
// \partial J / \partial y^i
- grad[globalDof][1] += weight * (-partA1 * sin(theta) + partA3 * cos(theta)) * shapeGrad[dof][0][0];
+ grad[globalDof][1] += weight * (-partA1 * sin(theta) + partA3 * cos(theta)) * shapeGrad[dof][0];
// \partial J / \partial \theta^i
- grad[globalDof][2] += weight * (B * theta_s * shapeGrad[dof][0][0]
+ grad[globalDof][2] += weight * (B * theta_s * shapeGrad[dof][0]
+ partA1 * (-x_s * sin(theta) - y_s * cos(theta)) * shapeFunction[dof]
+ partA3 * ( x_s * cos(theta) - y_s * sin(theta)) * shapeFunction[dof]);
@@ -358,13 +378,15 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
}
-template <class FunctionSpaceType, int polOrd>
-double Dune::RodAssembler<FunctionSpaceType, polOrd>::
+template <class GridType, int polOrd>
+double Dune::RodAssembler<GridType, polOrd>::
computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
{
- const int maxlevel = grid_->maxlevel();
+ const int maxlevel = grid_->maxLevel();
double energy = 0;
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(maxlevel);
+
if (sol.size()!=grid_->size(maxlevel, gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
@@ -375,13 +397,17 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
for (; it!=endIt; ++it) {
// Extract local solution on this element
- const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
- const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();
-
+ // const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
+// const int numOfBaseFct = baseSet.getNumberOfBaseFunctions();
+ const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->geometry().type(), elementOrder);
+ int numOfBaseFct = baseSet.size();
+
FieldVector<double, blocksize> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+ //localSolution[i] = sol[functionSpace_.mapToGlobal(*it,i)];
+ localSolution[i] = sol[indexSet.template subIndex<gridDim>(*it,i)];
// Get quadrature rule
const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().type(), polOrd);
@@ -391,7 +417,7 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
// Local position of the quadrature point
const FieldVector<double,gridDim>& quadPos = quad[pt].position();
- const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverse(quadPos);
+ const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = it->geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement;
@@ -399,17 +425,19 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
// ///////////////////////////////////////
// Compute deformation gradient
// ///////////////////////////////////////
- Array<JacobianRange> shapeGrad(numOfBaseFct);
+ Array<FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
for (int dof=0; dof<numOfBaseFct; dof++) {
- baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+ //baseSet.jacobian(dof, quadPos, shapeGrad[dof]);
+ for (int i=0; i<gridDim; i++)
+ shapeGrad[dof][i] = baseSet[dof].evaluateDerivative(0,i,quadPos);
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
// multiply with jacobian inverse
- JacobianRange tmp(0);
- inv.umv(shapeGrad[dof][0], tmp[0]);
- shapeGrad[dof][0] = tmp[0];
+ FieldVector<double,gridDim> tmp(0);
+ inv.umv(shapeGrad[dof], tmp);
+ shapeGrad[dof] = tmp;
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
@@ -417,17 +445,17 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
}
// Get the value of the shape functions
- RangeType shapeFunction[2];
+ double shapeFunction[2];
for(int i=0; i<2; i++)
- baseSet.eval(i,quadPos,shapeFunction[i]);
+ shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
// //////////////////////////////////
// Interpolate
// //////////////////////////////////
- double x_s = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*shapeGrad[1][0][0];
- double y_s = localSolution[0][1]*shapeGrad[0][0][0] + localSolution[1][1]*shapeGrad[1][0][0];
- double theta_s = localSolution[0][2]*shapeGrad[0][0][0] + localSolution[1][2]*shapeGrad[1][0][0];
+ double x_s = localSolution[0][0]*shapeGrad[0][0] + localSolution[1][0]*shapeGrad[1][0];
+ double y_s = localSolution[0][1]*shapeGrad[0][0] + localSolution[1][1]*shapeGrad[1][0];
+ double theta_s = localSolution[0][2]*shapeGrad[0][0] + localSolution[1][2]*shapeGrad[1][0];
double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1];
@@ -443,7 +471,6 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
+ A1 * partA1 * partA1
+ A3 * partA3 * partA3);
-// printf("energy %g\n", energy);
}
}
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 839f89e7..057e8bc8 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -11,39 +11,25 @@ namespace Dune
/** \brief The FEM operator for an extensible, shearable rod
*/
- template <class FunctionSpaceType, int polOrd>
+ template <class GridType, int polOrd>
class RodAssembler {
- //! The grid
- typedef typename FunctionSpaceType::GridType GridType;
-
typedef typename GridType::template Codim<0>::Entity EntityType;
typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
- typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
-
//! Dimension of the grid. This needs to be one!
enum { gridDim = GridType::dimension };
+ enum { elementOrder = 1};
+
//! Each block is x, y, theta
enum { blocksize = 3 };
//!
typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock;
- //! ???
- typedef typename FunctionSpaceType::JacobianRangeType JacobianRange;
-
- //! ???
- typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
- typedef typename FunctionSpaceType::RangeType RangeType;
-
- /** \todo Does actually belong into the base class */
const GridType* grid_;
- /** \todo Does actually belong into the base class */
- const FunctionSpaceType& functionSpace_;
-
/** \brief Material constants */
double B;
double A1;
@@ -52,10 +38,9 @@ namespace Dune
public:
//! ???
- RodAssembler(const FunctionSpaceType &f) :
- functionSpace_(f)
+ RodAssembler(const GridType &grid) :
+ grid_(&grid)
{
- grid_ = &f.grid();
B = 1;
A1 = 1;
A3 = 1;
diff --git a/staticrod.cc b/staticrod.cc
index c2d122dc..853ff7e6 100644
--- a/staticrod.cc
+++ b/staticrod.cc
@@ -1,20 +1,14 @@
#include <config.h>
+//#define DUNE_EXPRESSIONTEMPLATES
#include <dune/grid/onedgrid.hh>
-#include <dune/fem/lagrangebase.hh>
-#include <dune/grid/common/gridpart.hh>
-
#include <dune/istl/io.hh>
-//#include "../common/boundarytreatment.hh"
#include "../common/boundarypatch.hh"
#include <dune/common/bitfield.hh>
-//#include "../common/readbitfield.hh"
#include "src/rodassembler.hh"
-//#include "../common/linearipopt.hh"
-
#include "../common/projectedblockgsstep.hh"
#include "../contact/src/contactmmgstep.hh"
@@ -22,9 +16,10 @@
#include "../common/geomestimator.hh"
#include "../common/energynorm.hh"
-#include <dune/common/configparser.hh>
#include "src/rodwriter.hh"
+#include <dune/common/configparser.hh>
+
// Choose a solver
//#define IPOPT
//#define GAUSS_SEIDEL
@@ -92,7 +87,7 @@ int main (int argc, char *argv[]) try
// Problem settings
const int numRodBaseElements = parameterSet.get("numRodBaseElements", int(0));
-
+
// ///////////////////////////////////////
// Create the two grids
// ///////////////////////////////////////
@@ -103,7 +98,7 @@ int main (int argc, char *argv[]) try
for (int i=0; i<maxLevel; i++)
rod.globalRefine(1);
- int maxlevel = rod.maxlevel();
+ int maxlevel = rod.maxLevel();
int numRodElements = rod.size(maxlevel, 0);
@@ -119,23 +114,6 @@ int main (int argc, char *argv[]) try
}
}
- // //////////////////////////////////////////////////////////
- // Create discrete function spaces
- // //////////////////////////////////////////////////////////
-
- typedef FunctionSpace < double , double, 1, 1 > RodFuncSpace;
- typedef LevelGridPart<RodGridType> RodGridPartType;
- typedef LagrangeDiscreteFunctionSpace < RodFuncSpace, RodGridPartType, 1> RodFuncSpaceType;
-
- Array<RodGridPartType*> rodGridPart(maxlevel+1);
- Array<const RodFuncSpaceType*> rodFuncSpace(maxlevel+1);
-
- for (int i=0; i<maxlevel+1; i++) {
- rodGridPart[i] = new RodGridPartType(rod, i);
- rodFuncSpace[i] = new RodFuncSpaceType(*rodGridPart[i]);
- }
-
-
// ////////////////////////////////////////////////////////////
// Create solution and rhs vectors
// ////////////////////////////////////////////////////////////
@@ -145,17 +123,17 @@ int main (int argc, char *argv[]) try
VectorType corr;
MatrixType hessianMatrix;
- RodAssembler<RodFuncSpaceType, 4> rodAssembler(*rodFuncSpace[maxlevel]);
+ RodAssembler<RodGridType,4> rodAssembler(rod);
+
rodAssembler.setParameters(1, 100, 100);
MatrixIndexSet indices(numRodElements+1, numRodElements+1);
rodAssembler.getNeighborsPerVertex(indices);
indices.exportIdx(hessianMatrix);
-
- rhs.resize(rodFuncSpace[maxlevel]->size());
- x.resize(rodFuncSpace[maxlevel]->size());
- corr.resize(rodFuncSpace[maxlevel]->size());
+ rhs.resize(rod.size(maxlevel,1));
+ x.resize(rod.size(maxlevel,1));
+ corr.resize(rod.size(maxlevel,1));
// Initial solution
x = 0;
@@ -243,7 +221,7 @@ int main (int argc, char *argv[]) try
ProjectedBlockGSStep<MatrixType, VectorType> presmoother;
ProjectedBlockGSStep<MatrixType, VectorType> postsmoother;
- ContactMMGStep<MatrixType, VectorType, RodFuncSpaceType > contactMMGStep(maxlevel+1);
+ ContactMMGStep<MatrixType, VectorType, RodGridType > contactMMGStep(maxlevel+1);
contactMMGStep.setMGType(mu, nu1, nu2);
contactMMGStep.dirichletNodes_ = &dirichletNodes;
@@ -257,7 +235,8 @@ int main (int argc, char *argv[]) try
contactMMGStep.mgTransfer_.resize(maxlevel);
for (int i=0; i<contactMMGStep.mgTransfer_.size(); i++){
TruncatedMGTransfer<VectorType>* newTransferOp = new TruncatedMGTransfer<VectorType>;
- newTransferOp->setup(*rodFuncSpace[i], *rodFuncSpace[i+1]);
+ //newTransferOp->setup(*rodFuncSpace[i], *rodFuncSpace[i+1]);
+ newTransferOp->setup(rod,i,i+1);
contactMMGStep.mgTransfer_[i] = newTransferOp;
}
@@ -308,6 +287,7 @@ int main (int argc, char *argv[]) try
corr = 0;
//std::cout <<"Solution: " << x << std::endl;
+ //exit(0);
rodAssembler.assembleGradient(x, rhs);
rodAssembler.assembleMatrix(x, hessianMatrix);
@@ -407,11 +387,9 @@ int main (int argc, char *argv[]) try
//break;
} while (loadFactor < 1);
-
-
-
} catch (Exception e) {
std::cout << e << std::endl;
+
}
--
GitLab