Skip to content
Snippets Groups Projects
Commit 33fd5621 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Remove the RodAssembler class 

GeodesicFEAssembler can replace it without problems.
parent ac686e40
No related branches found
No related tags found
1 merge request!67Modernize rod code
Expand all Hide whitespace changes
Inline Side-by-side
dune/gfe/rodassembler.cc deleted 100644 → 0
+ 0
− 430
View file @ ac686e40
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/localfunctions/lagrange/p1.hh>
#include <dune/gfe/rodlocalstiffness.hh>
template <class Basis>
void RodAssembler<Basis,3>::
assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
using namespace Dune;
if (sol.size()!=this->basis_.size())
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
grad = 0;
// A view on the FE basis on a single element
auto localView = this->basis_.localView();
// Loop over all elements
for (const auto& element : Dune::elements(this->basis_.gridView()))
{
localView.bind(element);
// A 1d grid has two vertices
static const int nDofs = 2;
// Extract local solution
std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)];
// Assemble local gradient
std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
this->localStiffness_->assembleGradient(localView,
localSolution,
localGradient);
// Add to global gradient
for (int i=0; i<nDofs; i++)
grad[localView.index(i)] += localGradient[i];
}
}
template <class GridView>
void RodAssembler<GridView,2>::
assembleMatrix(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix)
{
Dune::MatrixIndexSet neighborsPerVertex;
this->getNeighborsPerVertex(neighborsPerVertex);
matrix = 0;
Dune::Matrix<MatrixBlock> mat;
for (const auto& element : Dune::elements(this->basis_.gridView()))
{
const int numOfBaseFct = 2;
// Extract local solution
std::vector<RigidBodyMotion<double,2> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[this->basis_.index(element,i)];
// setup matrix
getLocalMatrix(element, localSolution, mat);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
int row = this->basis_.index(element,i);
for (int j=0; j<numOfBaseFct; j++ ) {
int col = this->basis_.index(element,j);
matrix[row][col] += mat[i][j];
}
}
}
}
template <class GridView>
void RodAssembler<GridView,2>::
getLocalMatrix( EntityType &entity,
const std::vector<RigidBodyMotion<double,2> >& localSolution,
Dune::Matrix<MatrixBlock>& localMat) const
{
/* ndof is the number of vectors of the element */
int ndof = localSolution.size();
localMat.setSize(ndof,ndof);
localMat = 0;
Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
// Get quadrature rule
const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(entity.type(), 2);
/* Loop over all integration points */
for (int ip=0; ip<quad.size(); ip++) {
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[ip].position();
// calc Jacobian inverse before integration element is evaluated
const Dune::FieldMatrix<double,gridDim,gridDim>& inv = entity.geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = entity.geometry().integrationElement(quadPos);
/* Compute the weight of the current integration point */
double weight = quad[ip].weight() * integrationElement;
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
std::vector<Dune::FieldMatrix<double,1,gridDim> > referenceElementGradients(ndof);
localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
std::vector<Dune::FieldVector<double,gridDim> > shapeGrad(ndof);
// multiply with jacobian inverse
for (int dof=0; dof<ndof; dof++)
inv.mv(referenceElementGradients[dof][0], shapeGrad[dof]);
std::vector<Dune::FieldVector<double,1> > shapeFunction;
localFiniteElement.localBasis().evaluateFunction(quadPos,shapeFunction);
// //////////////////////////////////
// Interpolate
// //////////////////////////////////
double x_s = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
double y_s = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
double theta = localSolution[0].q.angle_*shapeFunction[0] + localSolution[1].q.angle_*shapeFunction[1];
for (int i=0; i<localMat.N(); i++) {
for (int j=0; j<localMat.M(); j++) {
// \partial J^2 / \partial x_i \partial x_j
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] * 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][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] * 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][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][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][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][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][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))
+ A3 * (x_s*cos(theta) - y_s*sin(theta)) * (x_s*cos(theta) - y_s*sin(theta))
- A3 * (x_s*sin(theta) + y_s*cos(theta) - 1) * (x_s*sin(theta) + y_s*cos(theta)));
}
}
}
#if 0
static int eleme = 0;
printf("********* Element %d **********\n", eleme++);
for (int row=0; row<matSize; row++) {
for (int rcomp=0; rcomp<dim; rcomp++) {
for (int col=0; col<matSize; col++) {
for (int ccomp=0; ccomp<dim; ccomp++)
std::cout << mat[row][col][rcomp][ccomp] << " ";
std::cout << " ";
}
std::cout << std::endl;
}
std::cout << std::endl;
}
exit(0);
#endif
}
template <class GridView>
void RodAssembler<GridView,2>::
assembleGradient(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
if (sol.size()!=this->basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
grad = 0;
// Loop over all elements
for (const auto& element : Dune::elements(this->basis_gridView()))
{
// Extract local solution on this element
Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
const int numOfBaseFct = localFiniteElement.localBasis().size();
RigidBodyMotion<double,2> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[this->basis_.index(element,i)];
// Get quadrature rule
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(element.type(), 2);
for (int pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
const Dune::FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = element.geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement;
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
std::vector<Dune::FieldMatrix<double,1,gridDim> > referenceElementGradients(numOfBaseFct);
localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
std::vector<Dune::FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
// multiply with jacobian inverse
for (int dof=0; dof<numOfBaseFct; dof++)
inv.mv(referenceElementGradients[dof][0], shapeGrad[dof]);
// Get the values of the shape functions
std::vector<Dune::FieldVector<double,1> > shapeFunction;
localFiniteElement.localBasis().evaluateFunction(quadPos,shapeFunction);
// //////////////////////////////////
// Interpolate
// //////////////////////////////////
double x_s = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
double y_s = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
double theta_s = localSolution[0].q.angle_*shapeGrad[0][0] + localSolution[1].q.angle_*shapeGrad[1][0];
double theta = localSolution[0].q.angle_*shapeFunction[0] + localSolution[1].q.angle_*shapeFunction[1];
// /////////////////////////////////////////////
// Sum it all up
// /////////////////////////////////////////////
double partA1 = A1 * (x_s * cos(theta) - y_s * sin(theta));
double partA3 = A3 * (x_s * sin(theta) + y_s * cos(theta) - 1);
for (int dof=0; dof<numOfBaseFct; dof++) {
int globalDof = this->basis_.index(element,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];
// \partial J / \partial y^i
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]
+ partA1 * (-x_s * sin(theta) - y_s * cos(theta)) * shapeFunction[dof]
+ partA3 * ( x_s * cos(theta) - y_s * sin(theta)) * shapeFunction[dof]);
}
}
}
}
template <class GridView>
double RodAssembler<GridView,2>::
computeEnergy(const std::vector<RigidBodyMotion<double,2> >& sol) const
{
double energy = 0;
if (sol.size()!=this->basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
// Loop over all elements
for (const auto& element : Dune::elements(this->basis_gridView()))
{
// Extract local solution on this element
Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement;
int numOfBaseFct = localFiniteElement.localBasis().size();
std::vector<RigidBodyMotion<double,2> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[this->basis_.index(element,i)];
// Get quadrature rule
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(element.type(), 2);
for (int pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
const Dune::FieldMatrix<double,1,1>& inv = element.geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = element.geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement;
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
std::vector<Dune::FieldMatrix<double,1,gridDim> > referenceElementGradients(numOfBaseFct);
localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceElementGradients);
std::vector<Dune::FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
// multiply with jacobian inverse
for (int dof=0; dof<numOfBaseFct; dof++)
inv.mv(referenceElementGradients[dof][0], shapeGrad[dof]);
// Get the value of the shape functions
std::vector<Dune::FieldVector<double,1> > shapeFunction;
localFiniteElement.localBasis().evaluateFunction(quadPos,shapeFunction);
// //////////////////////////////////
// Interpolate
// //////////////////////////////////
double x_s = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
double y_s = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
double theta_s = localSolution[0].q.angle_*shapeGrad[0][0] + localSolution[1].q.angle_*shapeGrad[1][0];
double theta = localSolution[0].q.angle_*shapeFunction[0] + localSolution[1].q.angle_*shapeFunction[1];
// /////////////////////////////////////////////
// Sum it all up
// /////////////////////////////////////////////
double partA1 = x_s * cos(theta) - y_s * sin(theta);
double partA3 = x_s * sin(theta) + y_s * cos(theta) - 1;
energy += 0.5 * weight * (B * theta_s * theta_s
+ A1 * partA1 * partA1
+ A3 * partA3 * partA3);
}
}
return energy;
}
dune/gfe/rodassembler.hh deleted 100644 → 0
+ 0
− 148
View file @ ac686e40
#ifndef DUNE_EXTENSIBLE_ROD_ASSEMBLER_HH
#define DUNE_EXTENSIBLE_ROD_ASSEMBLER_HH
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/gfe/localgeodesicfefdstiffness.hh>
#include "rigidbodymotion.hh"
#include "rodlocalstiffness.hh"
#include "geodesicfeassembler.hh"
/** \brief The FEM operator for an extensible, shearable rod in 3d
*/
template <class Basis, int spaceDim>
class RodAssembler
{
static_assert(spaceDim==2 || spaceDim==3,
"You can only instantiate the class RodAssembler for 2d and 3d spaces");
};
/** \brief The FEM operator for an extensible, shearable rod in 3d
*/
template <class Basis>
class RodAssembler<Basis,3> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >
{
typedef typename Basis::GridView GridView;
//! Dimension of the grid.
enum { gridDim = GridView::dimension };
static_assert(gridDim==1, "RodAssembler can only be used with one-dimensional grids!");
enum { elementOrder = 1};
//! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
enum { blocksize = 6 };
public:
//! ???
RodAssembler(const Basis& basis,
LocalGeodesicFEStiffness<Basis, RigidBodyMotion<double,3> >& localStiffness)
: GeodesicFEAssembler<Basis, RigidBodyMotion<double,3> >(basis,localStiffness)
{
std::vector<RigidBodyMotion<double,3> > referenceConfiguration(basis.size());
for (const auto vertex : Dune::vertices(basis.gridView()))
{
auto idx = basis.gridView().indexSet().index(vertex);
referenceConfiguration[idx].r[0] = 0;
referenceConfiguration[idx].r[1] = 0;
referenceConfiguration[idx].r[2] = vertex.geometry().corner(0)[0];
referenceConfiguration[idx].q = Rotation<double,3>::identity();
}
rodEnergy()->setReferenceConfiguration(referenceConfiguration);
}
auto rodEnergy()
{
// TODO: Does not work for other stiffness implementations
auto localFDStiffness = std::dynamic_pointer_cast<LocalGeodesicFEFDStiffness<Basis, RigidBodyMotion<double,3> > >(this->localStiffness_);
return const_cast<RodLocalStiffness<GridView,double>*>(dynamic_cast<const RodLocalStiffness<GridView,double>*>(localFDStiffness->localEnergy_));
}
std::vector<RigidBodyMotion<double,3> > getRefConfig()
{
return rodEnergy()->referenceConfiguration_;
}
virtual void assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const override;
}; // end class
/** \brief The FEM operator for a 2D extensible, shearable rod
*/
template <class Basis>
class RodAssembler<Basis,2> : public GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >
{
typedef typename Basis::GridView GridView;
typedef typename GridView::template Codim<0>::Entity EntityType;
//! Dimension of the grid.
enum { gridDim = GridView::dimension };
static_assert(gridDim==1, "RodAssembler can only be used with one-dimensional grids!");
enum { elementOrder = 1};
//! Each block is x, y, theta
enum { blocksize = 3 };
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
/** \brief Material constants */
double B;
double A1;
double A3;
public:
//! ???
RodAssembler(const GridView &gridView)
: GeodesicFEAssembler<Basis, RigidBodyMotion<double,2> >(gridView,nullptr)
{
B = 1;
A1 = 1;
A3 = 1;
}
~RodAssembler() {}
void setParameters(double b, double a1, double a3) {
B = b;
A1 = a1;
A3 = a3;
}
/** \brief Assemble the tangent stiffness matrix and the right hand side
*/
void assembleMatrix(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix);
virtual void assembleGradient(const std::vector<RigidBodyMotion<double,2> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const override;
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<RigidBodyMotion<double,2> >& sol) const override;
protected:
/** \brief Compute the element tangent stiffness matrix */
void getLocalMatrix( EntityType &entity,
const std::vector<RigidBodyMotion<double,2> >& localSolution,
Dune::Matrix<MatrixBlock>& mat) const;
}; // end class
#include "rodassembler.cc"
#endif
src/rod3d.cc
+ 27
− 9
View file @ 33fd5621
......@@ -13,9 +13,11 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/geodesicfeassembler.hh>
#include <dune/gfe/localgeodesicfefdstiffness.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/rodlocalstiffness.hh>
#include <dune/gfe/rotation.hh>
#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
typedef RigidBodyMotion<double,3> TargetSpace;
......@@ -118,16 +120,35 @@ int main (int argc, char *argv[]) try
dirichletNodes[0] = true;
dirichletNodes.back() = true;
//////////////////////////////////////////////
// Create the stress-free configuration
//////////////////////////////////////////////
auto localRodEnergy = std::make_shared<RodLocalStiffness<GridView, double> >(gridView,
A, J1, J2, E, nu);
std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(1));
for (const auto vertex : vertices(gridView))
{
auto idx = gridView.indexSet().index(vertex);
referenceConfiguration[idx].r[0] = 0;
referenceConfiguration[idx].r[1] = 0;
referenceConfiguration[idx].r[2] = vertex.geometry().corner(0)[0];
referenceConfiguration[idx].q = Rotation<double,3>::identity();
}
localRodEnergy->setReferenceConfiguration(referenceConfiguration);
// ///////////////////////////////////////////
// Create a solver for the rod problem
// ///////////////////////////////////////////
RodLocalStiffness<GridView,double> localStiffness(gridView,
A, J1, J2, E, nu);
LocalGeodesicFEFDStiffness<FEBasis,
TargetSpace> localStiffness(localRodEnergy.get());
LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localStiffness);
RodAssembler<FEBasis,3> rodAssembler(gridView, localFDStiffness);
GeodesicFEAssembler<FEBasis,TargetSpace> rodAssembler(gridView, localStiffness);
RiemannianTrustRegionSolver<FEBasis,RigidBodyMotion<double,3> > rodSolver;
......@@ -161,9 +182,6 @@ int main (int argc, char *argv[]) try
// //////////////////////////////
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "rod3d-result");
BlockVector<FieldVector<double, 6> > strain(x.size()-1);
rodAssembler.getStrain(x,strain);
} catch (Exception& e)
{
std::cout << e.what() << std::endl;
......
test/CMakeLists.txt
+ 0
− 1
View file @ 33fd5621
......@@ -9,7 +9,6 @@ set(TESTS
localprojectedfefunctiontest
orthogonalmatrixtest
rigidbodymotiontest
rodassemblertest
rotationtest
svdtest
targetspacetest
......
test/frameinvariancetest.cc
+ 22
− 11
View file @ 33fd5621
......@@ -5,14 +5,9 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/gfe/quaternion.hh>
#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/rodlocalstiffness.hh>
// Number of degrees of freedom:
// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
const int blocksize = 6;
using namespace Dune;
......@@ -74,14 +69,30 @@ int main (int argc, char *argv[]) try
rotatedX[i].q = rotation.mult(x[i].q);
}
RodLocalStiffness<GridView,double> localRodFirstOrderModel(gridView,
1,1,1,1e6,0.3);
RodLocalStiffness<GridView,double> localRodEnergy(gridView,
1,1,1,1e6,0.3);
std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(1));
for (const auto vertex : vertices(gridView))
{
auto idx = gridView.indexSet().index(vertex);
referenceConfiguration[idx].r[0] = 0;
referenceConfiguration[idx].r[1] = 0;
referenceConfiguration[idx].r[2] = vertex.geometry().corner(0)[0];
referenceConfiguration[idx].q = Rotation<double,3>::identity();
}
localRodEnergy.setReferenceConfiguration(referenceConfiguration);
LocalGeodesicFEFDStiffness<FEBasis,RigidBodyMotion<double,3> > localFDStiffness(&localRodFirstOrderModel);
auto localView = feBasis.localView();
localView.bind(*gridView.begin<0>());
RodAssembler<FEBasis,3> assembler(feBasis, localFDStiffness);
SolutionType localX = {x[0], x[1]};
SolutionType localRotatedX = {rotatedX[0], rotatedX[1]};
if (std::abs(assembler.computeEnergy(x) - assembler.computeEnergy(rotatedX)) > 1e-6)
if (std::abs(localRodEnergy.energy(localView, localX) - localRodEnergy.energy(localView, localRotatedX)) > 1e-6)
DUNE_THROW(Dune::Exception, "Rod energy not invariant under rigid body motions!");
} catch (Exception& e) {
......
test/rodassemblertest.cc deleted 100644 → 0
+ 0
− 587
View file @ ac686e40
This diff is collapsed.
test/rotationtest.cc
+ 0
− 127
View file @ 33fd5621
......@@ -5,12 +5,7 @@
#include <dune/common/fmatrix.hh>
#warning Do not include onedgrid.hh
#include <dune/grid/onedgrid.hh>
#include <dune/gfe/rotation.hh>
#warning Do not include rodlocalstiffness.hh
#include <dune/gfe/rodlocalstiffness.hh>
#include "valuefactory.hh"
using namespace Dune;
......@@ -94,123 +89,6 @@ void testDDExp()
}
}
void testDerivativeOfInterpolatedPosition()
{
std::array<Rotation<double,3>, 6> q;
FieldVector<double,3> xAxis = {1,0,0};
FieldVector<double,3> yAxis = {0,1,0};
FieldVector<double,3> zAxis = {0,0,1};
q[0] = Rotation<double,3>(xAxis, 0);
q[1] = Rotation<double,3>(xAxis, M_PI/2);
q[2] = Rotation<double,3>(yAxis, 0);
q[3] = Rotation<double,3>(yAxis, M_PI/2);
q[4] = Rotation<double,3>(zAxis, 0);
q[5] = Rotation<double,3>(zAxis, M_PI/2);
double eps = 1e-7;
for (int i=0; i<6; i++) {
for (int j=0; j<6; j++) {
for (int k=0; k<7; k++) {
double s = k/6.0;
std::array<Quaternion<double>,6> fdGrad;
// ///////////////////////////////////////////////////////////
// First: test the interpolated position
// ///////////////////////////////////////////////////////////
for (int l=0; l<3; l++) {
SkewMatrix<double,3> forward(FieldVector<double,3>(0));
SkewMatrix<double,3> backward(FieldVector<double,3>(0));
forward.axial()[l] += eps;
backward.axial()[l] -= eps;
fdGrad[l] = Rotation<double,3>::interpolate(q[i].mult(Rotation<double,3>::exp(forward)), q[j], s);
fdGrad[l] -= Rotation<double,3>::interpolate(q[i].mult(Rotation<double,3>::exp(backward)), q[j], s);
fdGrad[l] /= 2*eps;
fdGrad[3+l] = Rotation<double,3>::interpolate(q[i], q[j].mult(Rotation<double,3>::exp(forward)), s);
fdGrad[3+l] -= Rotation<double,3>::interpolate(q[i], q[j].mult(Rotation<double,3>::exp(backward)), s);
fdGrad[3+l] /= 2*eps;
}
// Compute analytical gradient
std::array<Quaternion<double>,6> grad;
RodLocalStiffness<OneDGrid,double>::interpolationDerivative(q[i], q[j], s, grad);
for (int l=0; l<6; l++) {
Quaternion<double> diff = fdGrad[l];
diff -= grad[l];
if (diff.two_norm() > 1e-6) {
std::cout << "Error in position " << l << ": fd: " << fdGrad[l]
<< " analytical: " << grad[l] << std::endl;
}
}
// ///////////////////////////////////////////////////////////
// Second: test the interpolated velocity vector
// ///////////////////////////////////////////////////////////
for (const double intervalLength : {1.0/3, 2.0/3, 3.0/3, 4.0/3, 5.0/3, 6.0/3, 7.0/3})
{
Dune::FieldVector<double,3> variation;
for (int m=0; m<3; m++) {
variation = 0;
variation[m] = eps;
fdGrad[m] = Rotation<double,3>::interpolateDerivative(q[i].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))),
q[j], s);
variation = 0;
variation[m] = -eps;
fdGrad[m] -= Rotation<double,3>::interpolateDerivative(q[i].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))),
q[j], s);
fdGrad[m] /= 2*eps;
variation = 0;
variation[m] = eps;
fdGrad[3+m] = Rotation<double,3>::interpolateDerivative(q[i], q[j].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))), s);
variation = 0;
variation[m] = -eps;
fdGrad[3+m] -= Rotation<double,3>::interpolateDerivative(q[i], q[j].mult(Rotation<double,3>::exp(SkewMatrix<double,3>(variation))), s);
fdGrad[3+m] /= 2*eps;
}
// Scale the finite difference gradient with the interval length
for (auto& parDer : fdGrad)
parDer /= intervalLength;
// Compute analytical velocity vector gradient
RodLocalStiffness<OneDGrid,double>::interpolationVelocityDerivative(q[i], q[j], s*intervalLength, intervalLength, grad);
for (int m=0; m<6; m++) {
Quaternion<double> diff = fdGrad[m];
diff -= grad[m];
if (diff.two_norm() > 1e-6) {
std::cout << "Error in velocity " << m
<< ": s = " << s << " of (" << intervalLength << ")"
<< " fd: " << fdGrad[m] << " analytical: " << grad[m] << std::endl;
}
}
}
}
}
}
}
void testRotation(Rotation<double,3> q)
{
// Make sure it really is a unit quaternion
......@@ -388,11 +266,6 @@ int main (int argc, char *argv[]) try
// //////////////////////////////////////////////
testDDExp();
// //////////////////////////////////////////////
// Test derivative of interpolated position
// //////////////////////////////////////////////
testDerivativeOfInterpolatedPosition();
return not passed;
} catch (Exception& e) {
......
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