From 9164962a81ddb83d5aa76b3906a844df19cc425c Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 29 May 2006 14:37:34 +0000
Subject: [PATCH] many fixes and a method to compute the strains
[[Imported from SVN: r817]]
---
src/rodassembler.cc | 227 ++++++++++++++++++++++++++++++++++++--------
src/rodassembler.hh | 3 +
2 files changed, 188 insertions(+), 42 deletions(-)
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index af4dc3f9..1f8584e5 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -3,7 +3,7 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/common/matrix.hh>
-#include <dune/quadrature/quadraturerules.hh>
+#include <dune/grid/common/quadraturerules.hh>
#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
@@ -93,13 +93,6 @@ assembleMatrix(const std::vector<Configuration>& sol,
}
-#if 0
- // temporary: make identity
- matrix = 0;
- for (int i=0; i<matrix.N(); i++)
- for (int j=0; j<6; j++)
- matrix[i][i][j][j] = 1;
-#endif
}
@@ -123,7 +116,6 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<matSize; j++)
localMat[i][j] = 0;
- //const BaseFunctionSetType & baseSet = this->functionSpace_.getBaseFunctionSet( entity );
const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
= Dune::LagrangeShapeFunctions<double, double, gridDim>::general(entity.geometry().type(), elementOrder);
@@ -200,6 +192,7 @@ getLocalMatrix( EntityType &entity,
}
Quaternion<double> dq_dvij_dvkl[2][3][2][3];
+ Quaternion<double> dq_dvij_dvkl_ds[2][3][2][3];
for (int i=0; i<2; i++) {
for (int j=0; j<3; j++) {
@@ -208,10 +201,14 @@ getLocalMatrix( EntityType &entity,
for (int l=0; l<3; l++) {
- for (int m=0; m<3; m++)
+ for (int m=0; m<3; m++) {
dq_dvij_dvkl[i][j][k][l][m] = 0;
+ dq_dvij_dvkl_ds[i][j][k][l][m] = 0;
+ }
dq_dvij_dvkl[i][j][k][l][3] = -0.25 * (j==l) * shapeFunction[i] * shapeFunction[k];
+ dq_dvij_dvkl_ds[i][j][k][l][3] = -0.25 * (j==l) * shapeGrad[i] * shapeGrad[k];
+
}
}
@@ -278,7 +275,7 @@ getLocalMatrix( EntityType &entity,
for (int l=0; l<3; l++) {
- FieldMatrix<double,3,3> A;
+ FieldMatrix<double,4,4> A;
for (int a=0; a<4; a++)
for (int b=0; b<4; b++)
A[a][b] = (dq_dvij[k][l].mult(hatq))[a] * (dq_dvij[i][j].mult(hatq))[b]
@@ -312,7 +309,29 @@ getLocalMatrix( EntityType &entity,
}
+ // Get the derivative of the rotation at the quadrature point by interpolating in $H$
+ Quaternion<double> hatq_s;
+ hatq_s[0] = localSolution[0].q[0]*shapeGrad[0] + localSolution[1].q[0]*shapeGrad[1];
+ hatq_s[1] = localSolution[0].q[1]*shapeGrad[0] + localSolution[1].q[1]*shapeGrad[1];
+ hatq_s[2] = localSolution[0].q[2]*shapeGrad[0] + localSolution[1].q[2]*shapeGrad[1];
+ hatq_s[3] = localSolution[0].q[3]*shapeGrad[0] + localSolution[1].q[3]*shapeGrad[1];
+ FieldVector<double,3> u; // The Darboux vector
+ u[0] = 2 * ( hatq[3]*hatq_s[0] + hatq[2]*hatq_s[1] - hatq[1]*hatq_s[2] - hatq[0]*hatq_s[3]);
+ u[1] = 2 * (-hatq[2]*hatq_s[0] + hatq[3]*hatq_s[1] + hatq[0]*hatq_s[2] - hatq[1]*hatq_s[3]);
+ u[2] = 2 * ( hatq[1]*hatq_s[0] - hatq[0]*hatq_s[1] + hatq[3]*hatq_s[2] - hatq[2]*hatq_s[3]);
+
+ // Contains \partial q / \partial v^i_j at v = 0
+ double dum_dvij[3][2][3];
+
+ for (int i=0; i<2; i++)
+ for (int j=0; j<3; j++) {
+
+ for (int m=0; m<3; m++)
+ dum_dvij[m][i][j] = B(m, dq_dvij[i][j].mult(hatq))*hatq_s + B(m,hatq)*(dq_dvij_ds[i][j].mult(hatq));
+
+ }
+
// ///////////////////////////////////
// Sum it all up
// ///////////////////////////////////
@@ -335,7 +354,7 @@ getLocalMatrix( EntityType &entity,
+ A3 * shapeGrad[i] * hatq.director(2)[j] * shapeGrad[k] * hatq.director(2)[l]);
// \partial W^2 \partial v^i_j \partial r^k_l
- localMat[i][k][j][l+3] = weight
+ localMat[i][k][j][l+3] += weight
* (A1 * shapeGrad[k]*hatq.director(0)[l]*(r_s*dd_dvij[0][i][j])
+ A1 * (r_s*hatq.director(0)) * shapeGrad[k] * dd_dvij[0][i][j][l]
+ A2 * shapeGrad[k]*hatq.director(1)[l]*(r_s*dd_dvij[1][i][j])
@@ -346,7 +365,7 @@ getLocalMatrix( EntityType &entity,
localMat[i][k][j+3][l] = localMat[i][k][j][l+3];
// \partial W^2 \partial v^i_j \partial v^k_l
- localMat[i][k][j+3][l+3] = weight
+ localMat[i][k][j+3][l+3] += weight
* (A1 * (r_s * dd_dvij[0][k][l]) * (r_s * dd_dvij[0][i][j])
+ A1 * (r_s * hatq.director(0)) * (r_s * dd_dvij_dvkl[0][i][j][k][l])
+ A2 * (r_s * dd_dvij[1][k][l]) * (r_s * dd_dvij[1][i][j])
@@ -357,9 +376,30 @@ getLocalMatrix( EntityType &entity,
// ////////////////////////////////////////////
// The rotational part
// ////////////////////////////////////////////
+ // Stupid: I want those as an array
+ double K[3] = {K1, K2, K3};
// \partial W^2 \partial v^i_j \partial v^k_l
// All other derivatives are zero
+ for (int m=0; m<3; m++) {
+
+ double sum = dum_dvij[m][k][l] * (B(m,dq_dvij[i][j].mult(hatq)) * hatq_s);
+
+ sum += dum_dvij[m][k][l] * (B(m,hatq)*(dq_dvij_ds[i][j].mult(hatq) + dq_dvij[i][j].mult(hatq_s)));
+
+ sum += u[m] * (B(m, dq_dvij_dvkl[i][j][k][l].mult(hatq)) * hatq_s);
+
+ sum += u[m] * (B(m, dq_dvij[i][j].mult(hatq)) * dq_dvij_ds[k][l].mult(hatq));
+
+ sum += u[m] * (B(m, dq_dvij[k][l].mult(hatq)) *
+ (dq_dvij_ds[i][j].mult(hatq) + dq_dvij[i][j].mult(hatq_s)));
+
+ sum += u[m] * (B(m, hatq) *
+ (dq_dvij_dvkl_ds[i][j][k][l].mult(hatq) + (dq_dvij_dvkl[i][j][k][l].mult(hatq_s))));
+
+ localMat[i][k][j+3][l+3] += 2*weight *K[m] * sum;
+
+ }
}
@@ -371,32 +411,6 @@ getLocalMatrix( EntityType &entity,
}
-#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 GridType>
@@ -482,7 +496,7 @@ assembleGradient(const std::vector<Configuration>& sol,
hatq[3] = localSolution[0].q[3]*shapeFunction[0] + localSolution[1].q[3]*shapeFunction[1];
hatq.normalize();
- // Get the derivative of the rotation at the quadrature point by interpolating in $H$ and normalizing
+ // Get the derivative of the rotation at the quadrature point by interpolating in $H$
Quaternion<double> hatq_s;
hatq_s[0] = localSolution[0].q[0]*shapeGrad[0] + localSolution[1].q[0]*shapeGrad[1];
hatq_s[1] = localSolution[0].q[1]*shapeGrad[0] + localSolution[1].q[1]*shapeGrad[1];
@@ -721,7 +735,7 @@ computeEnergy(const std::vector<Configuration>& sol) const
//std::cout << "strain : " << v << std::endl;
- energy += 0.5*A1*v[0]*v[0] + 0.5*A2*v[1]*v[1] + 0.5*A3*(v[2]-1)*(v[2]-1);
+ energy += weight * 0.5*A1*v[0]*v[0] + 0.5*A2*v[1]*v[1] + 0.5*A3*(v[2]-1)*(v[2]-1);
// Part II: the bending and twisting energy
@@ -732,7 +746,7 @@ computeEnergy(const std::vector<Configuration>& sol) const
//std::cout << "Darboux vector : " << u << std::endl;
- energy += 0.5 * (K1*u[0]*u[0] + K2*u[1]*u[1] + K3*u[2]*u[2]);
+ energy += weight * 0.5 * (K1*u[0]*u[0] + K2*u[1]*u[1] + K3*u[2]*u[2]);
}
@@ -741,3 +755,132 @@ computeEnergy(const std::vector<Configuration>& sol) const
return energy;
}
+
+
+template <class GridType>
+void Dune::RodAssembler<GridType>::
+getStrain(const std::vector<Configuration>& sol,
+ BlockVector<FieldVector<double, blocksize> >& strain) const
+{
+ const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet();
+
+ if (sol.size()!=indexSet.size(gridDim))
+ DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
+
+ // Strain defined on each element
+ strain.resize(indexSet.size(0));
+
+ ElementLeafIterator it = grid_->template leafbegin<0>();
+ ElementLeafIterator endIt = grid_->template leafend<0>();
+
+ // Loop over all elements
+ for (; it!=endIt; ++it) {
+
+ // Extract local solution on this element
+ const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->geometry().type(), elementOrder);
+ int numOfBaseFct = baseSet.size();
+
+ Configuration localSolution[numOfBaseFct];
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localSolution[i] = sol[indexSet.template subIndex<gridDim>(*it,i)];
+
+ // Get quadrature rule
+ const int polOrd = 2;
+ const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(it->geometry().type(), polOrd);
+
+ for (int pt=0; pt<quad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const FieldVector<double,gridDim>& quadPos = quad[pt].position();
+
+ const FieldMatrix<double,1,1>& inv = it->geometry().jacobianInverseTransposed(quadPos);
+ const double integrationElement = it->geometry().integrationElement(quadPos);
+
+ double weight = quad[pt].weight() * integrationElement;
+
+ // ///////////////////////////////////////
+ // Compute deformation gradient
+ // ///////////////////////////////////////
+ std::vector<FieldVector<double,gridDim> > shapeGrad(numOfBaseFct);
+
+ for (int dof=0; dof<numOfBaseFct; 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
+ FieldVector<double,gridDim> tmp(0);
+ inv.umv(shapeGrad[dof], tmp);
+ shapeGrad[dof] = tmp;
+ //std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
+
+ }
+
+ // Get the value of the shape functions
+ double shapeFunction[2];
+ for(int i=0; i<2; i++)
+ shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
+
+ // //////////////////////////////////
+ // Interpolate
+ // //////////////////////////////////
+
+ FieldVector<double,3> r_s;
+ r_s[0] = localSolution[0].r[0]*shapeGrad[0][0] + localSolution[1].r[0]*shapeGrad[1][0];
+ r_s[1] = localSolution[0].r[1]*shapeGrad[0][0] + localSolution[1].r[1]*shapeGrad[1][0];
+ r_s[2] = localSolution[0].r[2]*shapeGrad[0][0] + localSolution[1].r[2]*shapeGrad[1][0];
+
+ // Get the rotation at the quadrature point by interpolating in $H$ and normalizing
+ Quaternion<double> q;
+ q[0] = localSolution[0].q[0]*shapeFunction[0] + localSolution[1].q[0]*shapeFunction[1];
+ q[1] = localSolution[0].q[1]*shapeFunction[0] + localSolution[1].q[1]*shapeFunction[1];
+ q[2] = localSolution[0].q[2]*shapeFunction[0] + localSolution[1].q[2]*shapeFunction[1];
+ q[3] = localSolution[0].q[3]*shapeFunction[0] + localSolution[1].q[3]*shapeFunction[1];
+
+ // The interpolated quaternion is not a unit quaternion anymore. We simply normalize
+ q.normalize();
+
+ // Get the derivative of the rotation at the quadrature point by interpolating in $H$
+ Quaternion<double> q_s;
+ q_s[0] = localSolution[0].q[0]*shapeGrad[0][0] + localSolution[1].q[0]*shapeGrad[1][0];
+ q_s[1] = localSolution[0].q[1]*shapeGrad[0][0] + localSolution[1].q[1]*shapeGrad[1][0];
+ q_s[2] = localSolution[0].q[2]*shapeGrad[0][0] + localSolution[1].q[2]*shapeGrad[1][0];
+ q_s[3] = localSolution[0].q[3]*shapeGrad[0][0] + localSolution[1].q[3]*shapeGrad[1][0];
+
+ // /////////////////////////////////////////////
+ // Sum it all up
+ // /////////////////////////////////////////////
+
+ // Part I: the shearing and stretching strain
+ //std::cout << "tangent : " << r_s << std::endl;
+ FieldVector<double,3> v;
+ v[0] = r_s * q.director(0); // shear strain
+ v[1] = r_s * q.director(1); // shear strain
+ v[2] = r_s * q.director(2); // stretching strain
+
+ //std::cout << "strain : " << v << std::endl;
+
+ // Part II: the Darboux vector
+
+ FieldVector<double,3> u;
+ u[0] = 2 * ( q[3]*q_s[0] + q[2]*q_s[1] - q[1]*q_s[2] - q[0]*q_s[3]);
+ u[1] = 2 * (-q[2]*q_s[0] + q[3]*q_s[1] + q[0]*q_s[2] - q[1]*q_s[3]);
+ u[2] = 2 * ( q[1]*q_s[0] - q[0]*q_s[1] + q[3]*q_s[2] - q[2]*q_s[3]);
+
+ // Sum it all up
+ int elementIdx = indexSet.index(*it);
+ strain[elementIdx][0] += weight * v[0];
+ strain[elementIdx][1] += weight * v[1];
+ strain[elementIdx][2] += weight * v[2];
+ strain[elementIdx][3] += weight * u[0];
+ strain[elementIdx][4] += weight * u[1];
+ strain[elementIdx][5] += weight * u[2];
+
+ }
+
+ }
+
+}
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 173a62c6..ba2ac41e 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -71,6 +71,9 @@ namespace Dune
double computeEnergy(const std::vector<Configuration>& sol) const;
void getNeighborsPerVertex(MatrixIndexSet& nb) const;
+
+ void getStrain(const std::vector<Configuration>& sol,
+ BlockVector<FieldVector<double, blocksize> >& strain) const;
protected:
--
GitLab