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

bugfixes, and moved some stuff to the elasticityhelpers.hh file 

[[Imported from SVN: r9268]]
parent a5bcb27b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/neohookeassembler.cc
+ 59
− 162
View file @ 0eba4543
...@@ -7,93 +7,8 @@ ...@@ -7,93 +7,8 @@
#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh> #include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
#include <dune/quadrature/quadraturerules.hh> #include <dune/quadrature/quadraturerules.hh>
#include <dune/grid/common/referenceelements.hh> #include <dune/grid/common/referenceelements.hh>
static double delta(int i, int j) { return i==j? 1.0: 0.0; }
template <int dim>
double E_du(const Dune::FieldMatrix<double,dim,dim>& u, int i, int j, int p, int q)
{
double result = delta(i,p)*delta(j,q) + delta(j,p)*delta(i,q);
result += ( delta(i,q)*u[p][j] + u[p][i]*delta(j,q) );
return 0.5*result;
}
// template <int i, int j, int p, int q, int r, int s>
// static inline double e_dudu()
// {
// return 0.5 * Delta<p,r>::delta
// * ( Delta<i,q>::delta*Delta<j,s>::delta + Delta<i,s>::delta*Delta<j,q>::delta );
// }
static inline double e_dudu(int i, int j, int p, int q, int r, int s)
{
return 0.5 * delta(p,r)* ( delta(i,q)*delta(j,s) + delta(i,s)*delta(j,q) );
}
// The Kronecker delta function
template <int i, int j>
struct Delta {
static int const delta = 0;
};
template <int i> #include "elasticityhelpers.hh"
struct Delta<i,i> {
static int const delta = 1;
};
// \parder {det I + u}{ u_pq }
template <int p, int qq>
double det_du_tmpl(const Dune::FieldMatrix<double,3,3>& u) {
int const q = qq+3;
return (Delta<(p+1)%3,(q+1)%3>::delta + u[(p+1)%3][(q+1)%3])
* (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3])
- (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3])
* (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3]);
}
// \parder {det I + u}{ u_pq }
template <int p, int q>
double det_du_tmpl(const Dune::FieldMatrix<double,2,2>& u) {
return Delta<p,q>::delta * (1+u[(p+1)%2][(q+1)%2])
- Delta<p,(q+1)%2>::delta * u[(p+1)%2][(q+1)%2];
// p = 0, q = 0; (1+u[1][1])
// p = 0, q = 1; -u[1][0]
// p = 1, q = 0; -u[0][1]
// p = 1, q = 1; (1+u[0][0])
//
}
template <int p, int qq, int r, int s>
double det_dudu_tmpl(const Dune::FieldMatrix<double,3,3>& u) {
int const q = qq+3;
return
Delta<(p+1)%3,r>::delta*Delta<(q+1)%3,s>::delta * (Delta<(p+2)%3,(q+2)%3>::delta+u[(p+2)%3][(q+2)%3])
+ Delta<(p+2)%3,r>::delta*Delta<(q+2)%3,s>::delta * (Delta<(p+1)%3,(q+1)%3>::delta+u[(p+1)%3][(q+1)%3])
- Delta<(p+1)%3,r>::delta*Delta<(q-1)%3,s>::delta * (Delta<(p+2)%3,(q-2)%3>::delta+u[(p+2)%3][(q-2)%3])
- Delta<(p+2)%3,r>::delta*Delta<(q-2)%3,s>::delta * (Delta<(p+1)%3,(q-1)%3>::delta+u[(p+1)%3][(q-1)%3]);
}
// parder det (I + u) partial u_pq \partial u_rs
//template <int p, int q, int r, int s>
double det_dudu_tmpl(const Dune::FieldMatrix<double,2,2>& u, int p, int q, int r, int s) {
return delta(p,q)*delta(r,(p+1)%2)*delta(s,(q+1)%2)
- delta(p,(q+1)%2)*delta(r,(p+1)%2)*delta(s,(q+1)%2);
}
double det_val(const Dune::FieldMatrix<double,3,3>& u) {
return (1+u[0][0])*(1+u[1][1])*(1+u[2][2])
+ u[0][1]*u[1][2]*u[2][0]
+ u[0][2]*u[1][0]*u[2][1]
- u[2][0]*(1+u[1][1])*u[0][2]
- u[2][1]*u[1][2]*(1+u[0][0])
- (1+u[2][2])*u[1][0]*u[0][1];
}
double det_val(const Dune::FieldMatrix<double,2,2>& u) {
return (1+u[0][0])*(1+u[1][1]) - u[1][0]*u[0][1];
}
template <class GridType, int polOrd> template <class GridType, int polOrd>
void Dune::NeoHookeAssembler<GridType, polOrd>:: void Dune::NeoHookeAssembler<GridType, polOrd>::
...@@ -186,9 +101,13 @@ assembleMatrix(const BlockVector<FieldVector<double, dim> >& sol, ...@@ -186,9 +101,13 @@ assembleMatrix(const BlockVector<FieldVector<double, dim> >& sol,
} }
// Add elements of local rhs to global rhs // Add elements of local rhs to global rhs
for (int i=0; i<numOfBaseFct; i++) for (int i=0; i<numOfBaseFct; i++) {
rhs[indexSet.template subIndex<dim>(*it,i)] += localRhs[i]; rhs[indexSet.template subIndex<dim>(*it,i)] += localRhs[i];
// if (indexSet.template subIndex<dim>(*it,i) == 22 && rhs.size()==51) {
// std::cout << "rhs: " << localRhs[i] << " together " << rhs[22] << std::endl;
// std::cout << "element :\n" << localRhs << std::endl;
// }
}
} }
} }
...@@ -273,17 +192,6 @@ getLocalMatrix( EntityType &entity, ...@@ -273,17 +192,6 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<dim; j++) for (int j=0; j<dim; j++)
nablaPhi[dof*ncomp+comp][i][j] = (i==comp) ? shapeGrads[dof][j] : 0; nablaPhi[dof*ncomp+comp][i][j] = (i==comp) ? shapeGrads[dof][j] : 0;
#if 0
if (localSolution.two_norm() > 0.1) {
for (int i=0; i<ndof*ncomp; i++)
std::cout << "nablaPhi " << i << std::endl << nablaPhi[i] << std::endl;
}
#endif
// for (int i=0; i<ndof*ncomp; i++) {
// printf("Deformation gradient %d\n", i);
// std::cout << nablaPhi[i] << std::endl;
// }
/****************************************************/ /****************************************************/
// The deformation gradient of the deformation // The deformation gradient of the deformation
// in formulas: F(i,j) = $\partial \phi_i / \partial x_j$ // in formulas: F(i,j) = $\partial \phi_i / \partial x_j$
...@@ -301,28 +209,26 @@ getLocalMatrix( EntityType &entity, ...@@ -301,28 +209,26 @@ getLocalMatrix( EntityType &entity,
} }
//defGrad[i][i] += 1;
} }
//const double J = defGrad.determinant(); // J = det (I + \nabla U)
const double J = det_val(parDerU); const double J = Elasticity::det_val(parDerU);
// /////////////////////////////////////////////////////// // ///////////////////////////////////////////////////////
// The material moduli C_ijkl // The material moduli C_ijkl
// /////////////////////////////////////////////////////// // ///////////////////////////////////////////////////////
FieldMatrix<double,dim,dim> parJparU; FieldMatrix<double,dim,dim> parJparU;
parJparU[0][0] = det_du_tmpl<0,0>(parDerU); parJparU[0][0] = Elasticity::det_du_tmpl<0,0>(parDerU);
parJparU[0][1] = det_du_tmpl<0,1>(parDerU); parJparU[0][1] = Elasticity::det_du_tmpl<0,1>(parDerU);
parJparU[1][0] = det_du_tmpl<1,0>(parDerU); parJparU[1][0] = Elasticity::det_du_tmpl<1,0>(parDerU);
parJparU[1][1] = det_du_tmpl<1,1>(parDerU); parJparU[1][1] = Elasticity::det_du_tmpl<1,1>(parDerU);
if (dim==3) { if (dim==3) {
parJparU[0][2] = det_du_tmpl<0,2>(parDerU); parJparU[0][2] = Elasticity::det_du_tmpl<0,2>(parDerU);
parJparU[1][2] = det_du_tmpl<1,2>(parDerU); parJparU[1][2] = Elasticity::det_du_tmpl<1,2>(parDerU);
parJparU[2][0] = det_du_tmpl<2,0>(parDerU); parJparU[2][0] = Elasticity::det_du_tmpl<2,0>(parDerU);
parJparU[2][1] = det_du_tmpl<2,1>(parDerU); parJparU[2][1] = Elasticity::det_du_tmpl<2,1>(parDerU);
parJparU[2][2] = det_du_tmpl<2,2>(parDerU); parJparU[2][2] = Elasticity::det_du_tmpl<2,2>(parDerU);
} }
double C[dim][dim][dim][dim]; double C[dim][dim][dim][dim];
...@@ -335,16 +241,23 @@ getLocalMatrix( EntityType &entity, ...@@ -335,16 +241,23 @@ getLocalMatrix( EntityType &entity,
for (int l=0; l<dim; l++) { for (int l=0; l<dim; l++) {
C[i][j][k][l] = 0;
double tr_E_dudu = 0; double tr_E_dudu = 0;
for (int m=0; m<dim; m++) for (int m=0; m<dim; m++)
tr_E_dudu += e_dudu(m,m,i,j,k,l); tr_E_dudu += Elasticity::e_dudu(m,m,i,j,k,l);
C[i][j][k][l] += lambda_*J/2 * (parJparU[i][j]*parJparU[k][l] + J*det_dudu_tmpl(parDerU,i,j,k,l)); #define LAURSEN
C[i][j][k][l] += -(lambda_/2 + mu_)/J/J * (J*det_dudu_tmpl(parDerU,i,j,k,l) - parJparU[i][j]*parJparU[k][l]); // Laursen
#ifdef LAURSEN
C[i][j][k][l] = lambda_/2 * (parJparU[i][j]*parJparU[k][l] + J*Elasticity::det_dudu_tmpl(parDerU,i,j,k,l));
C[i][j][k][l] += -(lambda_/2 + mu_)/J/J * (J*Elasticity::det_dudu_tmpl(parDerU,i,j,k,l) - parJparU[i][j]*parJparU[k][l]);
C[i][j][k][l] += mu_ * tr_E_dudu; C[i][j][k][l] += mu_ * tr_E_dudu;
#else
// Fischer / Wriggers
C[i][j][k][l] = lambda_*(parJparU[i][j]*parJparU[k][l] + (J-1)*Elasticity::det_dudu_tmpl(parDerU,i,j,k,l));
//C[i][j][k][l] += -(2*mu_)/J/J * (J*det_dudu_tmpl(parDerU,i,j,k,l) - parJparU[i][j]*parJparU[k][l]);
C[i][j][k][l] += (2*mu_)*Gamma_xx(J)*det_dudu_tmpl(parDerU,i,j,k,l);
C[i][j][k][l] += mu_ * tr_E_dudu;
#endif
} }
} }
...@@ -362,7 +275,7 @@ getLocalMatrix( EntityType &entity, ...@@ -362,7 +275,7 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<dim; ++j) { for (int j=0; j<dim; ++j) {
trE_du[i][j] = 0; trE_du[i][j] = 0;
for (int k=0; k<dim; k++) for (int k=0; k<dim; k++)
trE_du[i][j] += E_du(parDerU,k,k,i,j); trE_du[i][j] += Elasticity::E_du(parDerU,k,k,i,j);
} }
FieldMatrix<double,dim,dim> fu(0); FieldMatrix<double,dim,dim> fu(0);
...@@ -371,32 +284,30 @@ getLocalMatrix( EntityType &entity, ...@@ -371,32 +284,30 @@ getLocalMatrix( EntityType &entity,
for (int j=0; j<dim; j++) { for (int j=0; j<dim; j++) {
// Laursen
#ifdef LAURSEN
fu[i][j] = J * lambda_ / 2 * parJparU[i][j]; fu[i][j] = J * lambda_ / 2 * parJparU[i][j];
fu[i][j] -= (lambda_/2 + mu_)/J * parJparU[i][j]; fu[i][j] -= (lambda_/2 + mu_)/J * parJparU[i][j];
fu[i][j] += mu_ * trE_du[i][j]; fu[i][j] += mu_ * trE_du[i][j];
#else
// Fischer
fu[i][j] = lambda_ * (J-1) * parJparU[i][j];
fu[i][j] += 2* mu_ *Gamma_x(J) * parJparU[i][j];
fu[i][j] += mu_ * trE_du[i][j];
#endif
} }
} }
//std::cout << "fu: " << std::endl << fu << std::endl;
// loop over all entries of the element stiffness matrix // loop over all entries of the element stiffness matrix
// and add the contribution of the current quadrature point // and add the contribution of the current quadrature point
for (int i=0; i<ndof*ncomp; i++) { for (int i=0; i<ndof*ncomp; i++) {
if (counter==0 || counter==7) {
std::cout << "nablaPhi " << i <<std::endl << nablaPhi[i];
//std::cout << "nablaPhi " << j <<std::endl << nablaPhi[j];
}
for (int j=0; j<ndof*ncomp; j++) { for (int j=0; j<ndof*ncomp; j++) {
// Compute \nabla \phi_i * C * \nabla \phi_j // Compute \nabla \phi_i * C * \nabla \phi_j
// if (counter==0 || counter==7) {
// std::cout << "nablaPhi " << i <<std::endl << nablaPhi[i];
// std::cout << "nablaPhi " << j <<std::endl << nablaPhi[j];
// }
// First tensor product: // First tensor product:
FieldMatrix<double,dim,dim> tmp(0); FieldMatrix<double,dim,dim> tmp(0);
for (int k=0; k<dim; k++) for (int k=0; k<dim; k++)
...@@ -408,8 +319,6 @@ getLocalMatrix( EntityType &entity, ...@@ -408,8 +319,6 @@ getLocalMatrix( EntityType &entity,
// Second inner tensor product. // Second inner tensor product.
for (int k=0; k<dim; k++) for (int k=0; k<dim; k++)
for (int l=0; l<dim; l++) { for (int l=0; l<dim; l++) {
// if (counter==0)
// printf("%d %d adding %g\n", i%dim, j%dim, nablaPhi[i][k][l] * tmp[k][l] * weight);
mat[i/dim][j/dim][i%dim][j%dim] += nablaPhi[i][k][l] * tmp[k][l] * weight; mat[i/dim][j/dim][i%dim][j%dim] += nablaPhi[i][k][l] * tmp[k][l] * weight;
} }
...@@ -525,8 +434,7 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const ...@@ -525,8 +434,7 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
for (int j=0; j<dim; j++) { for (int j=0; j<dim; j++) {
defGrad[i][j] = (i==j); defGrad[i][j] = 0;
//defGrad[i][j] = 0;
for (int k=0; k<numOfBaseFct; k++) for (int k=0; k<numOfBaseFct; k++)
defGrad[i][j] += localSolution[k][i]*shape_grads[k][j]; defGrad[i][j] += localSolution[k][i]*shape_grads[k][j];
...@@ -535,31 +443,13 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const ...@@ -535,31 +443,13 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
} }
// ///////////////////////////////////////
// Compute nonlinear strain tensor
// ///////////////////////////////////////
/** \todo Apparently we don't even need the full tensor but only its trace */
FieldMatrix<double,dim,dim> E;
// E = 1/2 (C - I) = 1/2 (F^T F - I)
for (int i=0; i<dim; i++) {
for (int j=0; j<dim; j++) {
E[i][j] = 0;
for (int k=0; k<dim; k++)
E[i][j] += defGrad[k][i] * defGrad[k][j];
E[i][j] -= (i==j);
E[i][j] /= 2;
}
}
// ///////////////////////////////////////////// // /////////////////////////////////////////////
// The trace of the nonlinear strain tensor // The trace of the nonlinear strain tensor
// ///////////////////////////////////////////// // /////////////////////////////////////////////
FieldMatrix<double,dim,dim> E;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
E[i][j] = Elasticity::E_val(defGrad, i, j);
double trE = 0; double trE = 0;
for (int i=0; i<dim; i++) for (int i=0; i<dim; i++)
...@@ -569,11 +459,18 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const ...@@ -569,11 +459,18 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
// Sum up result // Sum up result
// /////////////////////////////////////// // ///////////////////////////////////////
const double J = defGrad.determinant(); const double J = Elasticity::det_val(defGrad);//defGrad.determinant();
// printf("weight: %g, lambda: %g mu: %g J: %g, trE: %g, adding %g\n",
// quad[pt].weight(), lambda_, mu_, J, trE,
// weight * (lambda_ * (J*J-1)/4 - (lambda_/2 + mu_ ) * std::log(J) + mu_*trE));
// Laursen
#ifdef LAURSEN
energy += weight * (lambda_ * (J*J-1)/4 - (lambda_/2 + mu_ ) * std::log(J) + mu_*trE); energy += weight * (lambda_ * (J*J-1)/4 - (lambda_/2 + mu_ ) * std::log(J) + mu_*trE);
//energy += weight * (lambda_ * (J*-1)*(J-1)/4 - (lambda_/2+mu_) * std::log(J) + mu_*trE); #else
// Fischer
energy += weight * (lambda_ * (J-1)*(J-1)/2 + 2*mu_ * Gamma(J) + mu_*trE);
#endif
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment