diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 3ab8a6dba641eb72912312eedd99e5ad103bbae6..416266b48f113f301de06a9f2208e55fa6203063 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -51,14 +51,6 @@ void RodAssembler<GridType>::
assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix) const
{
- // ////////////////////////////////////////////////////
- // Create local assembler
- // ////////////////////////////////////////////////////
-
- Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
- Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
-
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
Dune::MatrixIndexSet neighborsPerVertex;
@@ -81,13 +73,8 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
- localStiffness.localReferenceConfiguration_.resize(numOfBaseFct);
-
- for (int i=0; i<numOfBaseFct; i++)
- localStiffness.localReferenceConfiguration_[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)];
-
// setup matrix
- localStiffness.assemble(*it, localSolution);
+ this->localStiffness_->assemble(*it, localSolution);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
@@ -97,7 +84,7 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
for (int j=0; j<numOfBaseFct; j++ ) {
int col = indexSet.subIndex(*it,j,gridDim);
- matrix[row][col] += localStiffness.mat(i,j);
+ matrix[row][col] += this->localStiffness_->mat(i,j);
}
}
@@ -119,14 +106,6 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
if (sol.size()!=grid_->size(maxlevel, gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
- // ////////////////////////////////////////////////////
- // Create local assembler
- // ////////////////////////////////////////////////////
-
- Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
- Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
-
grad.resize(sol.size());
grad = 0;
@@ -145,17 +124,10 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
- // Extract local reference configuration
- std::vector<RigidBodyMotion<3> > localReferenceConfiguration(nDofs);
-
- for (int i=0; i<nDofs; i++)
- localReferenceConfiguration[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)];
-
// Assemble local gradient
std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
- localStiffness.localReferenceConfiguration_ = localReferenceConfiguration;
- localStiffness.assembleGradient(*it, localSolution, localGradient);
+ this->localStiffness_->assembleGradient(*it, localSolution, localGradient);
// Add to global gradient
for (int i=0; i<nDofs; i++)
@@ -179,15 +151,6 @@ computeEnergy(const std::vector<RigidBodyMotion<3> >& sol) const
if (sol.size()!=indexSet.size(gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
- // ////////////////////////////////////////////////////
- // Create local assembler
- // ////////////////////////////////////////////////////
-
- Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
- Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
-
- std::vector<RigidBodyMotion<3> > localReferenceConfiguration(2);
std::vector<RigidBodyMotion<3> > localSolution(2);
ElementLeafIterator it = grid_->template leafbegin<0>();
@@ -196,15 +159,10 @@ computeEnergy(const std::vector<RigidBodyMotion<3> >& sol) const
// Loop over all elements
for (; it!=endIt; ++it) {
- for (int i=0; i<2; i++) {
-
- localReferenceConfiguration[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)];
+ for (int i=0; i<2; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
- }
-
- localStiffness.localReferenceConfiguration_ = localReferenceConfiguration;
- energy += localStiffness.energy(*it, localSolution);
+ energy += this->localStiffness_->energy(*it, localSolution);
}
@@ -225,14 +183,6 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
if (sol.size()!=indexSet.size(gridDim))
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
- // ////////////////////////////////////////////////////
- // Create local assembler
- // ////////////////////////////////////////////////////
-
- Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
- Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
-
// Strain defined on each element
strain.resize(indexSet.size(0));
strain = 0;
@@ -266,7 +216,7 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
double weight = quad[pt].weight() * it->geometry().integrationElement(quadPos);
- FieldVector<double,blocksize> localStrain = localStiffness.getStrain(localSolution, *it, quad[pt].position());
+ FieldVector<double,blocksize> localStrain = dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->getStrain(localSolution, *it, quad[pt].position());
// Sum it all up
strain[elementIdx].axpy(weight, localStrain);
@@ -295,13 +245,13 @@ getStress(const std::vector<RigidBodyMotion<3> >& sol,
// Get reference strain
Dune::BlockVector<Dune::FieldVector<double, blocksize> > referenceStrain;
- getStrain(referenceConfiguration_, referenceStrain);
+ getStrain(dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->referenceConfiguration_, referenceStrain);
// Linear diagonal constitutive law
for (size_t i=0; i<stress.size(); i++) {
for (int j=0; j<3; j++) {
- stress[i][j] = (stress[i][j] - referenceStrain[i][j]) * A_[j];
- stress[i][j+3] = (stress[i][j+3] - referenceStrain[i][j+3]) * K_[j];
+ stress[i][j] = (stress[i][j] - referenceStrain[i][j]) * dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->A_[j];
+ stress[i][j+3] = (stress[i][j+3] - referenceStrain[i][j+3]) * dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->K_[j];
}
}
}
@@ -347,29 +297,23 @@ getResultantForce(const BoundaryPatch<GridType>& boundary,
// Compute force across this boundary face
// //////////////////////////////////////////////
- // Create local assembler
-
- Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
- Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
- RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
-
double pos = nIt->intersectionSelfLocal().corner(0);
Dune::array<RigidBodyMotion<3>,2> localSolution = {sol[indexSet.template subIndex<1>(*eIt,0)],
- sol[indexSet.template subIndex<1>(*eIt,1)]};
- Dune::array<RigidBodyMotion<3>,2> localRefConf = {referenceConfiguration_[indexSet.template subIndex<1>(*eIt,0)],
- referenceConfiguration_[indexSet.template subIndex<1>(*eIt,1)]};
+ sol[indexSet.template subIndex<1>(*eIt,1)]};
+ Dune::array<RigidBodyMotion<3>,2> localRefConf = {dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->referenceConfiguration_[indexSet.template subIndex<1>(*eIt,0)],
+ dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->referenceConfiguration_[indexSet.template subIndex<1>(*eIt,1)]};
- FieldVector<double, blocksize> strain = localStiffness.getStrain(localSolution, *eIt, pos);
- FieldVector<double, blocksize> referenceStrain = localStiffness.getStrain(localRefConf, *eIt, pos);
+ FieldVector<double, blocksize> strain = dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->getStrain(localSolution, *eIt, pos);
+ FieldVector<double, blocksize> referenceStrain = dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->getStrain(localRefConf, *eIt, pos);
FieldVector<double,3> localStress;
for (int i=0; i<3; i++)
- localStress[i] = (strain[i] - referenceStrain[i]) * A_[i];
+ localStress[i] = (strain[i] - referenceStrain[i]) * dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->A_[i];
FieldVector<double,3> localTorque;
for (int i=0; i<3; i++)
- localTorque[i] = (strain[i+3] - referenceStrain[i+3]) * K_[i];
+ localTorque[i] = (strain[i+3] - referenceStrain[i+3]) * dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->K_[i];
// Transform stress given with respect to the basis given by the three directors to
// the canonical basis of R^3
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 407a52a78e40e5fc2290b449483cdefa1c1fd7a4..f20fa53412a29e74e2d64fdf2f82917c805a7010 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -38,13 +38,6 @@ class RodAssembler : public GeodesicFEAssembler<typename GridType::LeafGridView,
public:
const GridType* grid_;
- /** \brief Material constants */
- Dune::array<double,3> K_;
- Dune::array<double,3> A_;
-
- /** \brief The stress-free configuration */
- std::vector<RigidBodyMotion<3> > referenceConfiguration_;
-
public:
//! ???
@@ -54,11 +47,7 @@ class RodAssembler : public GeodesicFEAssembler<typename GridType::LeafGridView,
localStiffness),
grid_(&grid)
{
- // Set dummy material parameters
- K_[0] = K_[1] = K_[2] = 1;
- A_[0] = A_[1] = A_[2] = 1;
-
- referenceConfiguration_.resize(grid.size(gridDim));
+ std::vector<RigidBodyMotion<3> > referenceConfiguration(grid.size(gridDim));
typename GridType::template Codim<gridDim>::LeafIterator it = grid.template leafbegin<gridDim>();
typename GridType::template Codim<gridDim>::LeafIterator endIt = grid.template leafend<gridDim>();
@@ -67,45 +56,13 @@ class RodAssembler : public GeodesicFEAssembler<typename GridType::LeafGridView,
int idx = grid.leafIndexSet().index(*it);
- referenceConfiguration_[idx].r[0] = 0;
- referenceConfiguration_[idx].r[1] = 0;
- referenceConfiguration_[idx].r[2] = it->geometry().corner(0)[0];
- referenceConfiguration_[idx].q = Rotation<3,double>::identity();
+ referenceConfiguration[idx].r[0] = 0;
+ referenceConfiguration[idx].r[1] = 0;
+ referenceConfiguration[idx].r[2] = it->geometry().corner(0)[0];
+ referenceConfiguration[idx].q = Rotation<3,double>::identity();
}
- }
-
- void setParameters(double k1, double k2, double k3,
- double a1, double a2, double a3) {
- K_[0] = k1;
- K_[1] = k2;
- K_[2] = k3;
- A_[0] = a1;
- A_[1] = a2;
- A_[2] = a3;
- }
-
- /** \brief Set shape constants and material parameters
- \param A The rod section area
- \param J1, J2 The geometric moments (Flächenträgheitsmomente)
- \param E Young's modulus
- \param nu Poisson number
- */
- void setShapeAndMaterial(double A, double J1, double J2, double E, double nu)
- {
- // shear modulus
- double G = E/(2+2*nu);
-
- K_[0] = E * J1;
- K_[1] = E * J2;
- K_[2] = G * (J1 + J2);
-
- A_[0] = G * A;
- A_[1] = G * A;
- A_[2] = E * A;
-
- //printf("%g %g %g %g %g %g\n", K_[0], K_[1], K_[2], A_[0], A_[1], A_[2]);
- //exit(0);
+ dynamic_cast<RodLocalStiffness<typename GridType::LeafGridView, double>* >(this->localStiffness_)->setReferenceConfiguration(referenceConfiguration);
}
/** \brief Assemble the tangent stiffness matrix
diff --git a/src/rodlocalstiffness.hh b/src/rodlocalstiffness.hh
index d24c26639790786c082e08c1bf486dca2176eba2..7fd857629c5df926287a4a15dd530c080e690c07 100644
--- a/src/rodlocalstiffness.hh
+++ b/src/rodlocalstiffness.hh
@@ -30,7 +30,9 @@ class RodLocalStiffness
enum {bendingQuadOrder = 2};
public:
- std::vector<RigidBodyMotion<3> > localReferenceConfiguration_;
+
+ /** \brief The stress-free configuration */
+ std::vector<RigidBodyMotion<3> > referenceConfiguration_;
public:
@@ -54,12 +56,12 @@ public:
Dune::array<double,3> K_;
Dune::array<double,3> A_;
- //! Default Constructor
- RodLocalStiffness ()
- {}
+ GridView gridView_;
- //! Default Constructor
- RodLocalStiffness (const Dune::array<double,3>& K, const Dune::array<double,3>& A)
+ //! Constructor
+ RodLocalStiffness (const GridView& gridView,
+ const Dune::array<double,3>& K, const Dune::array<double,3>& A)
+ : gridView_(gridView)
{
for (int i=0; i<3; i++) {
K_[i] = K[i];
@@ -67,30 +69,64 @@ public:
}
}
- void assembleBoundaryCondition (const Entity& e, int k=1)
+ /** \brief Constructor setting shape constants and material parameters
+ \param A The rod section area
+ \param J1, J2 The geometric moments (Flächenträgheitsmomente)
+ \param E Young's modulus
+ \param nu Poisson number
+ */
+ RodLocalStiffness (const GridView& gridView,
+ double A, double J1, double J2, double E, double nu)
+ : gridView_(gridView)
{
- DUNE_THROW(Dune::NotImplemented, "!");
+ // shear modulus
+ double G = E/(2+2*nu);
+
+ K_[0] = E * J1;
+ K_[1] = E * J2;
+ K_[2] = G * (J1 + J2);
+
+ A_[0] = G * A;
+ A_[1] = G * A;
+ A_[2] = E * A;
}
+
+ void setReferenceConfiguration(const std::vector<RigidBodyMotion<3> >& referenceConfiguration) {
+ referenceConfiguration_ = referenceConfiguration;
+ }
+
virtual RT energy (const Entity& e,
const std::vector<RigidBodyMotion<3> >& localSolution) const;
+ /** \brief Assemble the element gradient of the energy functional */
+ void assembleGradient(const Entity& element,
+ const std::vector<RigidBodyMotion<3> >& solution,
+ std::vector<Dune::FieldVector<double,6> >& gradient) const;
+
+ Dune::FieldVector<double, 6> getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
+ const Entity& element,
+ const Dune::FieldVector<double,1>& pos) const;
+
+protected:
+
+ void getLocalReferenceConfiguration(const Entity& element,
+ std::vector<RigidBodyMotion<3> >& localReferenceConfiguration) const {
+
+ int numOfBaseFct = element.template count<dim>();
+ localReferenceConfiguration.resize(numOfBaseFct);
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localReferenceConfiguration[i] = referenceConfiguration_[gridView_.indexSet().subIndex(element,i,dim)];
+ }
+
static void interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
Dune::array<Quaternion<double>,6>& grad);
static void interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad);
- Dune::FieldVector<double, 6> getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
- const Entity& element,
- const Dune::FieldVector<double,1>& pos) const;
-
- /** \brief Assemble the element gradient of the energy functional */
- void assembleGradient(const Entity& element,
- const std::vector<RigidBodyMotion<3> >& solution,
- std::vector<Dune::FieldVector<double,6> >& gradient) const;
-
template <class T>
static Dune::FieldVector<T,3> darboux(const Rotation<3,T>& q, const Dune::FieldVector<T,4>& q_s)
{
@@ -109,11 +145,13 @@ template <class GridType, class RT>
RT RodLocalStiffness<GridType, RT>::
energy(const Entity& element,
const std::vector<RigidBodyMotion<3> >& localSolution
- //, const std::vector<RigidBodyMotion<3> >& localReferenceConfiguration,
) const
{
RT energy = 0;
+ std::vector<RigidBodyMotion<3> > localReferenceConfiguration;
+ getLocalReferenceConfiguration(element, localReferenceConfiguration);
+
// ///////////////////////////////////////////////////////////////////////////////
// The following two loops are a reduced integration scheme. We integrate
// the transverse shear and extensional energy with a first-order quadrature
@@ -135,7 +173,7 @@ energy(const Entity& element,
Dune::FieldVector<double,6> strain = getStrain(localSolution, element, quadPos);
// The reference strain
- Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration_, element, quadPos);
+ Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
for (int i=0; i<3; i++)
energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
@@ -156,7 +194,7 @@ energy(const Entity& element,
Dune::FieldVector<double,6> strain = getStrain(localSolution, element, quadPos);
// The reference strain
- Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration_, element, quadPos);
+ Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
// Part II: the bending and twisting energy
for (int i=0; i<3; i++)
@@ -462,6 +500,9 @@ assembleGradient(const Entity& element,
{
using namespace Dune;
+ std::vector<RigidBodyMotion<3> > localReferenceConfiguration;
+ getLocalReferenceConfiguration(element, localReferenceConfiguration);
+
// Extract local solution on this element
const Dune::LagrangeShapeFunctionSet<double, double, 1> & baseSet
= Dune::LagrangeShapeFunctions<double, double, 1>::general(element.type(), 1); // first order
@@ -524,7 +565,7 @@ assembleGradient(const Entity& element,
FieldVector<double,blocksize> strain = getStrain(solution, element, quadPos);
// The reference strain
- FieldVector<double,blocksize> referenceStrain = getStrain(localReferenceConfiguration_, element, quadPos);
+ FieldVector<double,blocksize> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
// dd_dvij[m][i][j] = \parder {(d_k)_i} {q}
@@ -600,7 +641,7 @@ assembleGradient(const Entity& element,
FieldVector<double,blocksize> strain = getStrain(solution, element, quadPos);
// The reference strain
- FieldVector<double,blocksize> referenceStrain = getStrain(localReferenceConfiguration_, element, quadPos);
+ FieldVector<double,blocksize> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
// First derivatives of the position
array<Quaternion<double>,6> dq_dwij;