Skip to content
Snippets Groups Projects
Commit 2fd33391 authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

using new base function implementations 

[[Imported from SVN: r597]]
parent 0051fbce
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/rodassembler.cc
+ 113
− 86
View file @ 2fd33391
...@@ -5,21 +5,22 @@ ...@@ -5,21 +5,22 @@
#include <dune/quadrature/quadraturerules.hh> #include <dune/quadrature/quadraturerules.hh>
#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
template <class GridType, int polOrd>
template <class FunctionSpaceType, int polOrd> void Dune::RodAssembler<GridType, polOrd>::
void Dune::RodAssembler<FunctionSpaceType, polOrd>::
getNeighborsPerVertex(MatrixIndexSet& nb) const getNeighborsPerVertex(MatrixIndexSet& nb) const
{ {
const int gridDim = GridType::dimension; const int gridDim = GridType::dimension;
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
int i, j; int i, j;
int n = grid_->size(grid_->maxlevel(), gridDim); int n = grid_->size(grid_->maxLevel(), gridDim);
nb.resize(n, n); nb.resize(n, n);
ElementIterator it = grid_->template lbegin<0>( grid_->maxlevel() ); ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
ElementIterator endit = grid_->template lend<0> ( grid_->maxlevel() ); ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
for (; it!=endit; ++it) { for (; it!=endit; ++it) {
...@@ -27,8 +28,10 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const ...@@ -27,8 +28,10 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const
for (j=0; j<it->template count<gridDim>(); j++) { for (j=0; j<it->template count<gridDim>(); j++) {
int iIdx = it->template subIndex<gridDim>(i); //int iIdx = it->template subIndex<gridDim>(i);
int jIdx = it->template subIndex<gridDim>(j); //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); nb.add(iIdx, jIdx);
...@@ -40,26 +43,29 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const ...@@ -40,26 +43,29 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const
} }
template <class FunctionSpaceType, int polOrd> template <class GridType, int polOrd>
void Dune::RodAssembler<FunctionSpaceType, polOrd>:: void Dune::RodAssembler<GridType, polOrd>::
assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol, assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
BCRSMatrix<MatrixBlock>& matrix) BCRSMatrix<MatrixBlock>& matrix)
{ {
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
MatrixIndexSet neighborsPerVertex; MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex); getNeighborsPerVertex(neighborsPerVertex);
matrix = 0; matrix = 0;
ElementIterator it = grid_->template lbegin<0>( grid_->maxlevel() ); ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
ElementIterator endit = grid_->template lend<0> ( grid_->maxlevel() ); ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
Matrix<MatrixBlock> mat; Matrix<MatrixBlock> mat;
for( ; it != endit; ++it ) { for( ; it != endit; ++it ) {
const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it ); //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);
const int numOfBaseFct = baseSet.size();
mat.resize(numOfBaseFct, numOfBaseFct); mat.resize(numOfBaseFct, numOfBaseFct);
...@@ -67,7 +73,8 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -67,7 +73,8 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
BlockVector<FieldVector<double, blocksize> > localSolution(numOfBaseFct); BlockVector<FieldVector<double, blocksize> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++) 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 // setup matrix
getLocalMatrix( *it, localSolution, numOfBaseFct, mat); getLocalMatrix( *it, localSolution, numOfBaseFct, mat);
...@@ -75,11 +82,13 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -75,11 +82,13 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
// Add element matrix to global stiffness matrix // Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) { 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++ ) { 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]; matrix[row][col] += mat[i][j];
} }
...@@ -94,14 +103,15 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -94,14 +103,15 @@ assembleMatrix(const BlockVector<FieldVector<double, blocksize> >& sol,
template <class FunctionSpaceType, int polOrd> template <class GridType, int polOrd>
template <class MatrixType> template <class MatrixType>
void Dune::RodAssembler<FunctionSpaceType, polOrd>:: void Dune::RodAssembler<GridType, polOrd>::
getLocalMatrix( EntityType &entity, getLocalMatrix( EntityType &entity,
const BlockVector<FieldVector<double, blocksize> >& localSolution, const BlockVector<FieldVector<double, blocksize> >& localSolution,
const int matSize, MatrixType& localMat) const 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 */ /* ndof is the number of vectors of the element */
int ndof = matSize; int ndof = matSize;
...@@ -109,9 +119,10 @@ getLocalMatrix( EntityType &entity, ...@@ -109,9 +119,10 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<matSize; j++) for (int j=0; j<matSize; j++)
localMat[i][j] = 0; 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 // Get quadrature rule
const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(entity.geometry().type(), polOrd); const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(entity.geometry().type(), polOrd);
...@@ -122,7 +133,7 @@ getLocalMatrix( EntityType &entity, ...@@ -122,7 +133,7 @@ getLocalMatrix( EntityType &entity,
const FieldVector<double,gridDim>& quadPos = quad[ip].position(); const FieldVector<double,gridDim>& quadPos = quad[ip].position();
// calc Jacobian inverse before integration element is evaluated // 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); const double integrationElement = entity.geometry().integrationElement(quadPos);
/* Compute the weight of the current integration point */ /* Compute the weight of the current integration point */
...@@ -131,31 +142,33 @@ getLocalMatrix( EntityType &entity, ...@@ -131,31 +142,33 @@ getLocalMatrix( EntityType &entity,
/**********************************************/ /**********************************************/
/* compute gradients of the shape functions */ /* compute gradients of the shape functions */
/**********************************************/ /**********************************************/
Array<JacobianRange> shapeGrad(ndof); Array<FieldVector<double,gridDim> > shapeGrad(ndof);
for (int dof=0; dof<ndof; dof++) { 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 // multiply with jacobian inverse
JacobianRange tmp(0); FieldVector<double,gridDim> tmp(0);
inv.umv(shapeGrad[dof][0], tmp[0]); inv.umv(shapeGrad[dof], tmp);
shapeGrad[dof][0] = tmp[0]; shapeGrad[dof] = tmp;
} }
RangeType shapeFunction[matSize]; double shapeFunction[matSize];
for(int i=0; i<matSize; i++) 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 // Interpolate
// ////////////////////////////////// // //////////////////////////////////
double x_s = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*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][0] + localSolution[1][1]*shapeGrad[1][0][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][0] + localSolution[1][2]*shapeGrad[1][0][0];
double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1]; double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1];
...@@ -164,51 +177,51 @@ getLocalMatrix( EntityType &entity, ...@@ -164,51 +177,51 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<matSize; j++) { for (int j=0; j<matSize; j++) {
// \partial J^2 / \partial x_i \partial x_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)); * (A1 * cos(theta) * cos(theta) + A3 * sin(theta) * sin(theta));
// \partial J^2 / \partial x_i \partial y_j // \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); * (-A1 + A3) * sin(theta)* cos(theta);
// \partial J^2 / \partial x_i \partial theta_j // \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*sin(theta) + y_s*cos(theta)) * cos(theta)
- A1* (x_s*cos(theta) - y_s*sin(theta)) * sin(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*cos(theta) - y_s*sin(theta)) * sin(theta)
+A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta)); +A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta));
// \partial J^2 / \partial y_i \partial x_j // \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)); * (-A1 * sin(theta)* cos(theta) + A3 * cos(theta) * sin(theta));
// \partial J^2 / \partial y_i \partial y_j // \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)); * (A1 * sin(theta)*sin(theta) + A3 * cos(theta)*cos(theta));
// \partial J^2 / \partial y_i \partial theta_j // \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 * sin(theta) + y_s * cos(theta)) * sin(theta)
-A1 * (x_s * cos(theta) - y_s * sin(theta)) * cos(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 * cos(theta) - y_s * sin(theta)) * cos(theta)
-A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta)); -A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta));
// \partial J^2 / \partial theta_i \partial x_j // \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*sin(theta) + y_s*cos(theta)) * cos(theta)
- A1* (x_s*cos(theta) - y_s*sin(theta)) * sin(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*cos(theta) - y_s*sin(theta)) * sin(theta)
+A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta)); +A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * cos(theta));
// \partial J^2 / \partial theta_i \partial y_j // \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 * sin(theta) + y_s * cos(theta)) * sin(theta)
-A1 * (x_s * cos(theta) - y_s * sin(theta)) * cos(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 * cos(theta) - y_s * sin(theta)) * cos(theta)
-A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta)); -A3 * (x_s * sin(theta) + y_s * cos(theta) - 1) * sin(theta));
// \partial J^2 / \partial theta_i \partial theta_j // \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] 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*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)) + A1 * (x_s*cos(theta) - y_s*sin(theta)) * (-x_s*cos(theta)+ y_s*sin(theta))
...@@ -252,12 +265,13 @@ getLocalMatrix( EntityType &entity, ...@@ -252,12 +265,13 @@ getLocalMatrix( EntityType &entity,
} }
template <class FunctionSpaceType, int polOrd> template <class GridType, int polOrd>
void Dune::RodAssembler<FunctionSpaceType, polOrd>:: void Dune::RodAssembler<GridType, polOrd>::
assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol, assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
BlockVector<FieldVector<double, blocksize> >& grad) const 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)) if (sol.size()!=grid_->size(maxlevel, gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!"); DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
...@@ -272,13 +286,16 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -272,13 +286,16 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
for (; it!=endIt; ++it) { for (; it!=endIt; ++it) {
// Extract local solution on this element // Extract local solution on this element
const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it ); //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);
const int numOfBaseFct = baseSet.size();
FieldVector<double, blocksize> localSolution[numOfBaseFct]; FieldVector<double, blocksize> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++) 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 // Get quadrature rule
const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().type(), polOrd); const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().type(), polOrd);
...@@ -288,7 +305,7 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -288,7 +305,7 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
// Local position of the quadrature point // Local position of the quadrature point
const FieldVector<double,gridDim>& quadPos = quad[pt].position(); 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); const double integrationElement = it->geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement; double weight = quad[pt].weight() * integrationElement;
...@@ -296,31 +313,33 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -296,31 +313,33 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
// /////////////////////////////////////// // ///////////////////////////////////////
// Compute deformation gradient // Compute deformation gradient
// /////////////////////////////////////// // ///////////////////////////////////////
Array<JacobianRange> shapeGrad(numOfBaseFct); Array<FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
for (int dof=0; dof<numOfBaseFct; dof++) { 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 // multiply with jacobian inverse
JacobianRange tmp(0); FieldVector<double,gridDim> tmp(0);
inv.umv(shapeGrad[dof][0], tmp[0]); inv.umv(shapeGrad[dof], tmp);
shapeGrad[dof][0] = tmp[0]; shapeGrad[dof] = tmp;
} }
// Get the value of the shape functions // Get the value of the shape functions
RangeType shapeFunction[2]; double shapeFunction[2];
for(int i=0; i<2; i++) for(int i=0; i<2; i++)
baseSet.eval(i,quadPos,shapeFunction[i]); shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
// ////////////////////////////////// // //////////////////////////////////
// Interpolate // Interpolate
// ////////////////////////////////// // //////////////////////////////////
double x_s = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*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][0] + localSolution[1][1]*shapeGrad[1][0][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][0] + localSolution[1][2]*shapeGrad[1][0][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]; double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1];
...@@ -333,18 +352,19 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol, ...@@ -333,18 +352,19 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
for (int dof=0; dof<numOfBaseFct; dof++) { 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); //printf("globalDof: %d partA1: %g partA3: %g\n", globalDof, partA1, partA3);
// \partial J / \partial x^i // \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 // \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 // \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] + partA1 * (-x_s * sin(theta) - y_s * cos(theta)) * shapeFunction[dof]
+ partA3 * ( x_s * cos(theta) - y_s * sin(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, ...@@ -358,13 +378,15 @@ assembleGradient(const BlockVector<FieldVector<double, blocksize> >& sol,
} }
template <class FunctionSpaceType, int polOrd> template <class GridType, int polOrd>
double Dune::RodAssembler<FunctionSpaceType, polOrd>:: double Dune::RodAssembler<GridType, polOrd>::
computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
{ {
const int maxlevel = grid_->maxlevel(); const int maxlevel = grid_->maxLevel();
double energy = 0; double energy = 0;
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(maxlevel);
if (sol.size()!=grid_->size(maxlevel, gridDim)) if (sol.size()!=grid_->size(maxlevel, gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!"); DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
...@@ -375,13 +397,17 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const ...@@ -375,13 +397,17 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
for (; it!=endIt; ++it) { for (; it!=endIt; ++it) {
// Extract local solution on this element // Extract local solution on this element
const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it ); // const BaseFunctionSetType & baseSet = functionSpace_.getBaseFunctionSet( *it );
const int numOfBaseFct = baseSet.getNumberOfBaseFunctions(); // 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]; FieldVector<double, blocksize> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++) 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 // Get quadrature rule
const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().type(), polOrd); 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 ...@@ -391,7 +417,7 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
// Local position of the quadrature point // Local position of the quadrature point
const FieldVector<double,gridDim>& quadPos = quad[pt].position(); 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); const double integrationElement = it->geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement; double weight = quad[pt].weight() * integrationElement;
...@@ -399,17 +425,19 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const ...@@ -399,17 +425,19 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
// /////////////////////////////////////// // ///////////////////////////////////////
// Compute deformation gradient // Compute deformation gradient
// /////////////////////////////////////// // ///////////////////////////////////////
Array<JacobianRange> shapeGrad(numOfBaseFct); Array<FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
for (int dof=0; dof<numOfBaseFct; dof++) { 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; //std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
// multiply with jacobian inverse // multiply with jacobian inverse
JacobianRange tmp(0); FieldVector<double,gridDim> tmp(0);
inv.umv(shapeGrad[dof][0], tmp[0]); inv.umv(shapeGrad[dof], tmp);
shapeGrad[dof][0] = tmp[0]; shapeGrad[dof] = tmp;
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl; //std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
...@@ -417,17 +445,17 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const ...@@ -417,17 +445,17 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
} }
// Get the value of the shape functions // Get the value of the shape functions
RangeType shapeFunction[2]; double shapeFunction[2];
for(int i=0; i<2; i++) for(int i=0; i<2; i++)
baseSet.eval(i,quadPos,shapeFunction[i]); shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
// ////////////////////////////////// // //////////////////////////////////
// Interpolate // Interpolate
// ////////////////////////////////// // //////////////////////////////////
double x_s = localSolution[0][0]*shapeGrad[0][0][0] + localSolution[1][0]*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][0] + localSolution[1][1]*shapeGrad[1][0][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][0] + localSolution[1][2]*shapeGrad[1][0][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]; double theta = localSolution[0][2]*shapeFunction[0] + localSolution[1][2]*shapeFunction[1];
...@@ -443,7 +471,6 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const ...@@ -443,7 +471,6 @@ computeEnergy(const BlockVector<FieldVector<double, blocksize> >& sol) const
+ A1 * partA1 * partA1 + A1 * partA1 * partA1
+ A3 * partA3 * partA3); + A3 * partA3 * partA3);
// printf("energy %g\n", energy);
} }
} }
......
src/rodassembler.hh
+ 5
− 20
View file @ 2fd33391
...@@ -11,39 +11,25 @@ namespace Dune ...@@ -11,39 +11,25 @@ namespace Dune
/** \brief The FEM operator for an extensible, shearable rod /** \brief The FEM operator for an extensible, shearable rod
*/ */
template <class FunctionSpaceType, int polOrd> template <class GridType, int polOrd>
class RodAssembler { class RodAssembler {
//! The grid
typedef typename FunctionSpaceType::GridType GridType;
typedef typename GridType::template Codim<0>::Entity EntityType; typedef typename GridType::template Codim<0>::Entity EntityType;
typedef typename GridType::template Codim<0>::LevelIterator ElementIterator; typedef typename GridType::template Codim<0>::LevelIterator ElementIterator;
typedef typename FunctionSpaceType::BaseFunctionSetType BaseFunctionSetType;
//! Dimension of the grid. This needs to be one! //! Dimension of the grid. This needs to be one!
enum { gridDim = GridType::dimension }; enum { gridDim = GridType::dimension };
enum { elementOrder = 1};
//! Each block is x, y, theta //! Each block is x, y, theta
enum { blocksize = 3 }; enum { blocksize = 3 };
//! //!
typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock; 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_; const GridType* grid_;
/** \todo Does actually belong into the base class */
const FunctionSpaceType& functionSpace_;
/** \brief Material constants */ /** \brief Material constants */
double B; double B;
double A1; double A1;
...@@ -52,10 +38,9 @@ namespace Dune ...@@ -52,10 +38,9 @@ namespace Dune
public: public:
//! ??? //! ???
RodAssembler(const FunctionSpaceType &f) : RodAssembler(const GridType &grid) :
functionSpace_(f) grid_(&grid)
{ {
grid_ = &f.grid();
B = 1; B = 1;
A1 = 1; A1 = 1;
A3 = 1; A3 = 1;
......
staticrod.cc
+ 15
− 37
View file @ 2fd33391
#include <config.h> #include <config.h>
//#define DUNE_EXPRESSIONTEMPLATES
#include <dune/grid/onedgrid.hh> #include <dune/grid/onedgrid.hh>
#include <dune/fem/lagrangebase.hh>
#include <dune/grid/common/gridpart.hh>
#include <dune/istl/io.hh> #include <dune/istl/io.hh>
//#include "../common/boundarytreatment.hh"
#include "../common/boundarypatch.hh" #include "../common/boundarypatch.hh"
#include <dune/common/bitfield.hh> #include <dune/common/bitfield.hh>
//#include "../common/readbitfield.hh"
#include "src/rodassembler.hh" #include "src/rodassembler.hh"
//#include "../common/linearipopt.hh"
#include "../common/projectedblockgsstep.hh" #include "../common/projectedblockgsstep.hh"
#include "../contact/src/contactmmgstep.hh" #include "../contact/src/contactmmgstep.hh"
...@@ -22,9 +16,10 @@ ...@@ -22,9 +16,10 @@
#include "../common/geomestimator.hh" #include "../common/geomestimator.hh"
#include "../common/energynorm.hh" #include "../common/energynorm.hh"
#include <dune/common/configparser.hh>
#include "src/rodwriter.hh" #include "src/rodwriter.hh"
#include <dune/common/configparser.hh>
// Choose a solver // Choose a solver
//#define IPOPT //#define IPOPT
//#define GAUSS_SEIDEL //#define GAUSS_SEIDEL
...@@ -92,7 +87,7 @@ int main (int argc, char *argv[]) try ...@@ -92,7 +87,7 @@ int main (int argc, char *argv[]) try
// Problem settings // Problem settings
const int numRodBaseElements = parameterSet.get("numRodBaseElements", int(0)); const int numRodBaseElements = parameterSet.get("numRodBaseElements", int(0));
// /////////////////////////////////////// // ///////////////////////////////////////
// Create the two grids // Create the two grids
// /////////////////////////////////////// // ///////////////////////////////////////
...@@ -103,7 +98,7 @@ int main (int argc, char *argv[]) try ...@@ -103,7 +98,7 @@ int main (int argc, char *argv[]) try
for (int i=0; i<maxLevel; i++) for (int i=0; i<maxLevel; i++)
rod.globalRefine(1); rod.globalRefine(1);
int maxlevel = rod.maxlevel(); int maxlevel = rod.maxLevel();
int numRodElements = rod.size(maxlevel, 0); int numRodElements = rod.size(maxlevel, 0);
...@@ -119,23 +114,6 @@ int main (int argc, char *argv[]) try ...@@ -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 // Create solution and rhs vectors
// //////////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////////
...@@ -145,17 +123,17 @@ int main (int argc, char *argv[]) try ...@@ -145,17 +123,17 @@ int main (int argc, char *argv[]) try
VectorType corr; VectorType corr;
MatrixType hessianMatrix; MatrixType hessianMatrix;
RodAssembler<RodFuncSpaceType, 4> rodAssembler(*rodFuncSpace[maxlevel]); RodAssembler<RodGridType,4> rodAssembler(rod);
rodAssembler.setParameters(1, 100, 100); rodAssembler.setParameters(1, 100, 100);
MatrixIndexSet indices(numRodElements+1, numRodElements+1); MatrixIndexSet indices(numRodElements+1, numRodElements+1);
rodAssembler.getNeighborsPerVertex(indices); rodAssembler.getNeighborsPerVertex(indices);
indices.exportIdx(hessianMatrix); indices.exportIdx(hessianMatrix);
rhs.resize(rod.size(maxlevel,1));
rhs.resize(rodFuncSpace[maxlevel]->size()); x.resize(rod.size(maxlevel,1));
x.resize(rodFuncSpace[maxlevel]->size()); corr.resize(rod.size(maxlevel,1));
corr.resize(rodFuncSpace[maxlevel]->size());
// Initial solution // Initial solution
x = 0; x = 0;
...@@ -243,7 +221,7 @@ int main (int argc, char *argv[]) try ...@@ -243,7 +221,7 @@ int main (int argc, char *argv[]) try
ProjectedBlockGSStep<MatrixType, VectorType> presmoother; ProjectedBlockGSStep<MatrixType, VectorType> presmoother;
ProjectedBlockGSStep<MatrixType, VectorType> postsmoother; ProjectedBlockGSStep<MatrixType, VectorType> postsmoother;
ContactMMGStep<MatrixType, VectorType, RodFuncSpaceType > contactMMGStep(maxlevel+1); ContactMMGStep<MatrixType, VectorType, RodGridType > contactMMGStep(maxlevel+1);
contactMMGStep.setMGType(mu, nu1, nu2); contactMMGStep.setMGType(mu, nu1, nu2);
contactMMGStep.dirichletNodes_ = &dirichletNodes; contactMMGStep.dirichletNodes_ = &dirichletNodes;
...@@ -257,7 +235,8 @@ int main (int argc, char *argv[]) try ...@@ -257,7 +235,8 @@ int main (int argc, char *argv[]) try
contactMMGStep.mgTransfer_.resize(maxlevel); contactMMGStep.mgTransfer_.resize(maxlevel);
for (int i=0; i<contactMMGStep.mgTransfer_.size(); i++){ for (int i=0; i<contactMMGStep.mgTransfer_.size(); i++){
TruncatedMGTransfer<VectorType>* newTransferOp = new TruncatedMGTransfer<VectorType>; 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; contactMMGStep.mgTransfer_[i] = newTransferOp;
} }
...@@ -308,6 +287,7 @@ int main (int argc, char *argv[]) try ...@@ -308,6 +287,7 @@ int main (int argc, char *argv[]) try
corr = 0; corr = 0;
//std::cout <<"Solution: " << x << std::endl; //std::cout <<"Solution: " << x << std::endl;
//exit(0);
rodAssembler.assembleGradient(x, rhs); rodAssembler.assembleGradient(x, rhs);
rodAssembler.assembleMatrix(x, hessianMatrix); rodAssembler.assembleMatrix(x, hessianMatrix);
...@@ -407,11 +387,9 @@ int main (int argc, char *argv[]) try ...@@ -407,11 +387,9 @@ int main (int argc, char *argv[]) try
//break; //break;
} while (loadFactor < 1); } while (loadFactor < 1);
} catch (Exception e) { } catch (Exception e) {
std::cout << e << std::endl; std::cout << e << std::endl;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment