Skip to content
Snippets Groups Projects
Commit b4e0b0da authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

bugfix: the explicit expression for the Darboux vector given by Dichmann is... 

bugfix: the explicit expression for the Darboux vector given by Dichmann is already in local coordinates

[[Imported from SVN: r1105]]
parent dd2172f6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/rodassembler.cc
+ 19
− 34
View file @ b4e0b0da
...@@ -307,9 +307,6 @@ getLocalMatrix( EntityPointer &entity, ...@@ -307,9 +307,6 @@ getLocalMatrix( EntityPointer &entity,
for (int i=0; i<4; i++) for (int i=0; i<4; i++)
hatq_s[i] = localSolution[0].q[i]*shapeGrad[0] + localSolution[1].q[i]*shapeGrad[1]; hatq_s[i] = localSolution[0].q[i]*shapeGrad[0] + localSolution[1].q[i]*shapeGrad[1];
// The Darboux vector
FieldVector<double,3> darbouxCan = darbouxCanonical(hatq, hatq_s);
// The current strain // The current strain
FieldVector<double,blocksize> strain = getStrain(globalSolution, entity, quadPos); FieldVector<double,blocksize> strain = getStrain(globalSolution, entity, quadPos);
...@@ -317,25 +314,22 @@ getLocalMatrix( EntityPointer &entity, ...@@ -317,25 +314,22 @@ getLocalMatrix( EntityPointer &entity,
FieldVector<double,blocksize> referenceStrain = getStrain(referenceConfiguration_, entity, quadPos); FieldVector<double,blocksize> referenceStrain = getStrain(referenceConfiguration_, entity, quadPos);
// Contains \partial u (canonical) / \partial v^i_j at v = 0 // Contains \partial u (canonical) / \partial v^i_j at v = 0
FieldVector<double,3> duCan_dvij[2][3]; FieldVector<double,3> du_dvij[2][3];
FieldVector<double,3> duLocal_dvij[2][3]; FieldVector<double,3> du_dvij_dvkl[2][3][2][3];
FieldVector<double,3> duCan_dvij_dvkl[2][3][2][3];
for (int i=0; i<2; i++) { for (int i=0; i<2; i++) {
for (int j=0; j<3; j++) { for (int j=0; j<3; j++) {
for (int m=0; m<3; m++) { for (int m=0; m<3; m++) {
duCan_dvij[i][j][m] = 2 * ( (hatq.mult(dq_dvj[j])).B(m)*hatq_s ) * shapeFunction[i]; //du_dvij[i][j][m] = 2 * (hatq.mult(dq_dvj[j])).B(m)*hatq_s * shapeFunction[i];
du_dvij[i][j][m] = (hatq.mult(dq_dvj[j])).B(m)*hatq_s;
du_dvij[i][j][m] *= 2 * shapeFunction[i];
Quaternion<double> tmp = dq_dvj[j]; Quaternion<double> tmp = dq_dvj[j];
tmp *= shapeFunction[i]; tmp *= shapeFunction[i];
duCan_dvij[i][j][m] += 2 * ( hatq.B(m)*(hatq.mult(dq_dvij_ds[i][j]) + hatq_s.mult(tmp))); du_dvij[i][j][m] += 2 * ( hatq.B(m)*(hatq.mult(dq_dvij_ds[i][j]) + hatq_s.mult(tmp)));
} }
// Don't fuse the two loops, we need the complete duCan_dvij[i][j]
for (int m=0; m<3; m++)
duLocal_dvij[i][j][m] = duCan_dvij[i][j] * hatq.director(m) + darbouxCan*dd_dvj[m][j]*shapeFunction[i];
for (int k=0; k<2; k++) { for (int k=0; k<2; k++) {
for (int l=0; l<3; l++) { for (int l=0; l<3; l++) {
...@@ -351,7 +345,7 @@ getLocalMatrix( EntityPointer &entity, ...@@ -351,7 +345,7 @@ getLocalMatrix( EntityPointer &entity,
Quaternion<double> tmp_ijkl = dq_dvj_dvl[j][l]; Quaternion<double> tmp_ijkl = dq_dvj_dvl[j][l];
tmp_ijkl *= shapeFunction[i] * shapeFunction[k]; tmp_ijkl *= shapeFunction[i] * shapeFunction[k];
duCan_dvij_dvkl[i][j][k][l] = 2 * ( (hatq.mult(tmp_ijkl)).B(m) * hatq_s) du_dvij_dvkl[i][j][k][l][m] = 2 * ( (hatq.mult(tmp_ijkl)).B(m) * hatq_s)
+ 2 * ( (hatq.mult(tmp_ij)).B(m) * (hatq.mult(dq_dvij_ds[k][l]) + hatq_s.mult(tmp_kl))) + 2 * ( (hatq.mult(tmp_ij)).B(m) * (hatq.mult(dq_dvij_ds[k][l]) + hatq_s.mult(tmp_kl)))
+ 2 * ( (hatq.mult(tmp_kl)).B(m) * (hatq.mult(dq_dvij_ds[i][j]) + hatq_s.mult(tmp_ij))) + 2 * ( (hatq.mult(tmp_kl)).B(m) * (hatq.mult(dq_dvij_ds[i][j]) + hatq_s.mult(tmp_ij)))
+ 2 * ( hatq.B(m) * (hatq.mult(dq_dvij_dvkl_ds[i][j][k][l]) + hatq_s.mult(tmp_ijkl))); + 2 * ( hatq.B(m) * (hatq.mult(dq_dvij_dvkl_ds[i][j][k][l]) + hatq_s.mult(tmp_ijkl)));
...@@ -409,12 +403,9 @@ getLocalMatrix( EntityPointer &entity, ...@@ -409,12 +403,9 @@ getLocalMatrix( EntityPointer &entity,
// All other derivatives are zero // All other derivatives are zero
//double sum = duLocal_dvij[k][l][m] * (duCan_dvij[i][j] * hatq.director(m) + darbouxCan*dd_dvj[m][j]*shapeFunction[i]); //double sum = duLocal_dvij[k][l][m] * (duCan_dvij[i][j] * hatq.director(m) + darbouxCan*dd_dvj[m][j]*shapeFunction[i]);
double sum = duLocal_dvij[k][l][m] * duLocal_dvij[i][j][m]; double sum = du_dvij[k][l][m] * du_dvij[i][j][m];
sum += (strain[m+3] - referenceStrain[m+3]) * (duCan_dvij_dvkl[i][j][k][l] * hatq.director(m) sum += (strain[m+3] - referenceStrain[m+3]) * du_dvij_dvkl[i][j][k][l][m];
+ duCan_dvij[i][j] * dd_dvj[m][l] * shapeFunction[k]
+ duCan_dvij[k][l] * dd_dvj[m][j] * shapeFunction[i]
+ darbouxCan * dd_dvij_dvkl[m][j][l] * shapeFunction[i] * shapeFunction[k]);
localMat[i][k][j+3][l+3] += weight *K_[m] * sum; localMat[i][k][j+3][l+3] += weight *K_[m] * sum;
...@@ -582,23 +573,17 @@ assembleGradient(const std::vector<Configuration>& sol, ...@@ -582,23 +573,17 @@ assembleGradient(const std::vector<Configuration>& sol,
double du_dvij_m; double du_dvij_m;
for (int t=0; t<3; t++) { du_dvij_m = (hatq.mult(dq_dvj[j])).B(m) * hatq_s;
du_dvij_m *= shapeFunction[i];
du_dvij_m = (hatq.mult(dq_dvj[j])).B(t) * hatq_s; Quaternion<double> tmp = dq_dvj[j];
du_dvij_m *= shapeFunction[i] * hatq.director(m)[t]; tmp *= shapeFunction[i];
Quaternion<double> tmp_ds = dq_dvj_ds[j];
Quaternion<double> tmp = dq_dvj[j]; tmp_ds *= shapeGrad[i];
tmp *= shapeFunction[i];
Quaternion<double> tmp_ds = dq_dvj_ds[j]; du_dvij_m += hatq.B(m) * (hatq.mult(tmp_ds) + hatq_s.mult(tmp));
tmp_ds *= shapeGrad[i];
du_dvij_m += hatq.B(t) * (hatq.mult(tmp_ds) + hatq_s.mult(tmp)) * hatq.director(m)[t];
du_dvij_m += hatq.B(t) * hatq_s * dd_dvj[m][j][t] * shapeFunction[i];
du_dvij_m *= 2;
} du_dvij_m *= 2;
grad[globalDof][3+j] += weight * K_[m] grad[globalDof][3+j] += weight * K_[m]
* (strain[m+3]-referenceStrain[m+3]) * du_dvij_m; * (strain[m+3]-referenceStrain[m+3]) * du_dvij_m;
......
src/rodassembler.hh
+ 5
− 19
View file @ b4e0b0da
...@@ -137,30 +137,16 @@ namespace Dune ...@@ -137,30 +137,16 @@ namespace Dune
const std::vector<Configuration>& globalSolution, const std::vector<Configuration>& globalSolution,
const int matSize, MatrixType& mat) const; const int matSize, MatrixType& mat) const;
template <class T> template <class T>
static FieldVector<T,3> darboux(const Quaternion<T>& q, const FieldVector<T,4>& q_s) static FieldVector<T,3> darboux(const Quaternion<T>& q, const FieldVector<T,4>& q_s)
{ {
FieldVector<double,3> uCanonical = darbouxCanonical(q, q_s); // The Darboux vector FieldVector<double,3> u; // The Darboux vector
FieldVector<double,3> u; u[0] = 2 * (q.B(0) * q_s);
u[0] = uCanonical*q.director(0); u[1] = 2 * (q.B(1) * q_s);
u[1] = uCanonical*q.director(1); u[2] = 2 * (q.B(2) * q_s);
u[2] = uCanonical*q.director(2);
return u;
}
template <class T> return u;
static FieldVector<T,3> darbouxCanonical(const Quaternion<T>& q, const FieldVector<T,4>& q_s)
{
FieldVector<double,3> uCanonical; // The Darboux vector
uCanonical[0] = 2 * (q.B(0) * q_s);
uCanonical[1] = 2 * (q.B(1) * q_s);
uCanonical[2] = 2 * (q.B(2) * q_s);
return uCanonical;
} }
static void getFirstDerivativesOfDirectors(const Quaternion<double>& q, static void getFirstDerivativesOfDirectors(const Quaternion<double>& q,
......
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