From d446c6f0bc495152628ab7869a86b4d433eec2a7 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Mon, 26 Jul 2021 23:37:12 +0200
Subject: [PATCH] Update L2-Error
---
src/dune-microstructure.cc | 420 ++++++++++++++++++++++++-------------
1 file changed, 269 insertions(+), 151 deletions(-)
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 0359ad4c..2510c3a2 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1378,6 +1378,151 @@ double evalSymGrad(const Basis& basis,
+
+template<class Basis, class Vector, class MatrixFunction>
+double computeL2ErrorTest(const Basis& basis,
+ Vector& coeffVector,
+ const double gamma,
+ const MatrixFunction& matrixFieldFunc
+ )
+{
+ auto error = 0.0;
+ constexpr int dim = 3;
+ constexpr int nCompo = 3;
+
+ auto localView = basis.localView();
+
+ auto matrixFieldGVF = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
+ auto matrixField = localFunction(matrixFieldGVF);
+
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
+
+ for (const auto& element : elements(basis.gridView()))
+ {
+ localView.bind(element);
+ matrixField.bind(element);
+ auto geometry = element.geometry();
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
+ const auto nSf = localFiniteElement.localBasis().size();
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+
+ for (const auto& quadPoint : quad)
+ {
+
+ const auto& quadPos = quadPoint.position();
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
+ const auto integrationElement = geometry.integrationElement(quadPoint.position());
+
+ std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
+ referenceGradients);
+
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+ // std::vector< VectorRT> gradients(referenceGradients.size());
+
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+
+ MatrixRT defGradientU(0), defGradientV(0);
+
+
+ for (size_t k=0; k < nCompo; k++)
+ for (size_t i=0; i < nSf; i++)
+
+ {
+ for (size_t l=0; l< nCompo; l++)
+ for (size_t j=0; j < nSf; j++ )
+ {
+
+ size_t localIdx1 = localView.tree().child(k).localIndex(i);
+ size_t globalIdx1 = localView.index(localIdx1);
+ size_t localIdx2 = localView.tree().child(l).localIndex(j);
+ size_t globalIdx2 = localView.index(localIdx2);
+
+ // (scaled) Deformation gradient of the ansatz basis function
+ defGradientU[k][0] = gradients[i][0]; // Y
+ defGradientU[k][1] = gradients[i][1]; //X2
+ defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
+
+ defGradientV[l][0] = gradients[j][0]; // Y
+ defGradientV[l][1] = gradients[j][1]; //X2
+ defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
+
+ // symmetric Gradient (Elastic Strains)
+ MatrixRT strainU(0), strainV(0);
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ {
+ strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
+ strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+ }
+
+ // Local energy density: stress times strain
+ double tmp = 0.0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
+
+
+ error += tmp * quadPoint.weight() * integrationElement;
+ }
+ }
+
+ for (size_t k=0; k < nCompo; k++)
+ for (size_t i=0; i < nSf; i++)
+ {
+ size_t localIdx1 = localView.tree().child(k).localIndex(i);
+ size_t globalIdx1 = localView.index(localIdx1);
+
+ // (scaled) Deformation gradient of the ansatz basis function
+ defGradientU[k][0] = gradients[i][0]; // Y
+ defGradientU[k][1] = gradients[i][1]; //X2
+ defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
+
+ // symmetric Gradient (Elastic Strains)
+ MatrixRT strainU(0);
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ {
+ strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
+ }
+
+ double tmp = 0.0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
+
+ error += tmp * quadPoint.weight() * integrationElement;
+ }
+
+ double tmp = 0.0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ tmp += matrixField(quadPos)[ii][jj] * matrixField(quadPos)[ii][jj];
+
+ error += tmp * quadPoint.weight() * integrationElement;
+
+ }
+ }
+ return sqrt(error);
+}
+
+
+
+
+
+
+
+
+
+
template<class Basis, class Vector, class MatrixFunction>
double computeL2Error(const Basis& basis,
Vector& coeffVector,
@@ -1400,9 +1545,6 @@ double computeL2Error(const Basis& basis,
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
-
-
for (const auto& element : elements(basis.gridView()))
{
@@ -1413,17 +1555,13 @@ double computeL2Error(const Basis& basis,
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
-
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
-
-
+
for (const auto& quadPoint : quad)
{
-
const auto& quadPos = quadPoint.position();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
const auto integrationElement = geometry.integrationElement(quadPoint.position());
@@ -1454,18 +1592,15 @@ double computeL2Error(const Basis& basis,
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-
-
-
+
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
+// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient //TEST
}
-
-
// Local energy density: stress times strain
double tmp = 0.0;
@@ -1473,16 +1608,12 @@ double computeL2Error(const Basis& basis,
for (int jj=0; jj<nCompo; jj++)
tmp+= std::pow(strainU[ii][jj] - matrixField(quadPos)[ii][jj],2);
// tmp+= std::pow((coeffVector[globalIdx]*strainU[ii][jj]) - matrixField(quadPos)[ii][jj] ,2);
-
-
-
+
error += tmp * quadPoint.weight() * integrationElement;
}
}
}
-
return sqrt(error);
-
}
@@ -1500,17 +1631,18 @@ int main(int argc, char *argv[])
ParameterTree parameterSet;
- // if (argc < 2)
- // ParameterTreeParser::readINITree("../src/cellsolver.parset", parameterSet);
- // ParameterTreeParser::readINITree(argv[1], parameterSet);
- // ParameterTreeParser::readOptions(argc, argv, parameterSet);
+// if (argc < 2)
+ ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
+
+// ParameterTreeParser::readINITree(argv[1], parameterSet);
+// ParameterTreeParser::readOptions(argc, argv, parameterSet);
// output setter
- // -- not set up --
+ std::string outputPath = parameterSet.get("outputPath", "../../outputs/output.txt");
std::fstream log;
std::fstream debugLog;
// log.open("output2.txt", std::ios::out);
- log.open("../../outputs/output.txt",std::ios::out);
+ log.open(outputPath ,std::ios::out);
debugLog.open("../../outputs/debugLog.txt",std::ios::out);
// log << "Writing this to a file. \n";
// log.close();
@@ -1525,14 +1657,7 @@ int main(int argc, char *argv[])
///////////////////////////////////
double gamma = parameterSet.get<double>("gamma",2000.0); // ratio dimension reduction to homogenization
- // Material parameter laminat
- // double E1 = parameterSet.get<double>("E1", 17e6); //material1
- // double nu1 = parameterSet.get<double>("nu1", 0.0);
- // double nu1 = parameterSet.get<double>("nu1", 0.3);
- // double E2 = parameterSet.get<double>("E2", 35e6); //material2
- // double nu2 = parameterSet.get<double>("nu2", 0.5);
-
- // // Plate geometry settings
+ // Plate geometry settings
double width = parameterSet.get<double>("width", 1.0); //geometry cell, cross section
// double len = parameterSet.get<double>("len", 10.0); //length
// double height = parameterSet.get<double>("h", 0.1); //rod thickness
@@ -1540,63 +1665,51 @@ int main(int argc, char *argv[])
///////////////////////////////////
- // ---Get Prestrain ---
+ // Get Prestrain Parameters
///////////////////////////////////
- unsigned int prestraintype = parameterSet.get<unsigned int>("imp", 1);
- double p1 = parameterSet.get<double>("prestress_pressure1", 1.0);
- double p2 = parameterSet.get<double>("prestress_pressure2", 2.0);
+ unsigned int prestraintype = parameterSet.get<unsigned int>("imp", 2);
+ double rho1 = parameterSet.get<double>("rho1", 1.0);
+ double alpha = parameterSet.get<double>("alpha", 2.0);
double theta = parameterSet.get<double>("theta",0.3);
- // double theta = 0.5;
- p1 = 1.0;
- double alpha = 2;
- p2 = alpha*1.0;
-
- prestraintype = 2;
- // prestraintype = 1;
+ double rho2 = alpha*1.0;
- auto prestrainImp = PrestrainImp(p1, p2, theta, width);
+ auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
auto B_Term = prestrainImp.getPrestrain(prestraintype);
- log << "prestrain imp: " << prestraintype << "\np1 = " << p1 << "\np2 = " << p2 << std::endl;
+ log << "prestrain imp: " << prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2 << std::endl;
log << "alpha: " << alpha << std::endl;
log << "gamma: " << gamma << std::endl;
-
- ///////////////////////////////////
- // --- Choose Solver ---
- // 1 : CG-Solver
- // 2 : GMRES
- // 3 : QR
- ///////////////////////////////////
- unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
- bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", true);
- bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", true);
- bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", true);
-
+
///////////////////////////////////
// Generate the grid
///////////////////////////////////
- using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
-
- // Domain : [0,1)^2 x (-1/2, 1/2)
- // FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
- // FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
-
+
+
+ //Corrector Problem Domain:
+ unsigned int cellDomain = parameterSet.get<unsigned int>("cellDomain", 1);
// (shifted) Domain (-1/2,1/2)^3
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
- int nE = 10;
- log << "Number of Elements: " << nE << std::endl;
- debugLog << "Number of Elements: " << nE << std::endl;
+ if (cellDomain == 2)
+ {
+ // Domain : [0,1)^2 x (-1/2, 1/2)
+ FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
+ FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
+ }
+
+// int nE = 10;
+// log << "Number of Elements: " << nE << std::endl;
+// debugLog << "Number of Elements: " << nE << std::endl;
- std::array<int,dim> nElements={nE,nE,nE};
+// std::array<int,dim> nElements={nE,nE,nE};
+ std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
// std::array<int,dim> nElements={100,100,1}; // TEST
-// printvector(std::cout, coeffT_1 , "coeffT_1", "--" );
-// printvector(std::cout, x_1 , "x_1", "--" );; //#Elements in each direction
-
- CellGridType grid_CE(lower,upper,nElements); //Corrector Problem Domain:
+ log << "Number of Elements in each direction: " << nElements << std::endl;
+ using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+ CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
@@ -1616,30 +1729,17 @@ int main(int argc, char *argv[])
using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
-
-
-
-// double beta = 1; //homogeneous case
- double beta = 2;
-
-
- double mu1 = 10;
-// double mu1 = 0.5*17e6;
-
- // double mu1 = 1000;
+ ///////////////////////////////////
+ // Get Material Parameters
+ ///////////////////////////////////
+ double beta = parameterSet.get<double>("beta",2.0);
+ double mu1 = parameterSet.get<double>("mu1",10.0);;
double mu2 = beta*mu1;
- // //
-
- // double mu1 = 0.5 * 1.0/(1.0 + nu1) * E1;
- // double lambda1 = nu1/(1.0-2.0*nu1) * 1.0/(1.0+nu1) * E1;
- double lambda1 = 0;
-
-
- // Create Lambda-Functions for material Parameters
- // auto muTerm = [mu1] (const Domain& z) {
- // return mu1;
- // };
-
+ double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+ double lambda2 = beta*lambda1;
+ ///////////////////////////////////
+ // Create Lambda-Functions for material Parameters depending on microstructure
+ ///////////////////////////////////
auto muTerm = [mu1, mu2, theta] (const Domain& z) {
if (abs(z[0]) >= (theta/2))
return mu1;
@@ -1650,24 +1750,39 @@ int main(int argc, char *argv[])
// return mu1;
};
-
-
-
-
- auto lambdaTerm = [lambda1] (const Domain& z) {
+ auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
+ if (abs(z[0]) >= (theta/2))
return lambda1;
+ else
+ return lambda2;
};
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE); // TODO needed here?!
auto muLocal = localFunction(muGridF);
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
auto lambdaLocal = localFunction(lambdaGridF);
log << "beta: " << beta << std::endl;
log << "material parameters: " << std::endl;
- log << " mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
- log << "\n lambda: " << lambda1 << std::endl;
+ log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+ log << " lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
+
+
+
+
+
+ ///////////////////////////////////
+ // --- Choose Solver ---
+ // 1 : CG-Solver
+ // 2 : GMRES
+ // 3 : QR
+ ///////////////////////////////////
+ unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+ bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
+ bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
+ bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
+
@@ -1692,9 +1807,9 @@ int main(int argc, char *argv[])
equivCounter++;
}
}
- std::cout << "equivCounter:" << equivCounter << std::endl;
- std::cout << "Number of Elements in each direction: " << nE << std::endl;
- std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
+// std::cout << "equivCounter:" << equivCounter << std::endl;
+// std::cout << "Number of Elements in each direction: " << nE << std::endl;
+// std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
auto perTest = periodicIndices.indexPairSet();
@@ -1707,15 +1822,14 @@ int main(int argc, char *argv[])
// blockedInterleaved())); // siehe Periodicbasistest.. funktioniert auch?
// flatInterleaved()));
-
- std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
- std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
-
+// std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
+// std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
+ log << "size feBasis: " << Basis_CE.size() << std::endl;
+
const int phiOffset = Basis_CE.size();
debugLog << "phiOffset: "<< phiOffset << std::endl;
-
// std::cout << Basis_CE.preBasis_.children << std::endl; // does not work
// std::cout << Basis_CE::preBasis_::children << std::endl;
// TEST
@@ -1731,7 +1845,7 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
- // Stiffness matrix and right hand side vector
+ // Data structures: Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
VectorCT load_alpha1, load_alpha2, load_alpha3;
MatrixCT stiffnessMatrix_CE;
@@ -1754,10 +1868,6 @@ int main(int argc, char *argv[])
debugLog << " substract_integralMean: " << substract_integralMean << std::endl;
-
- // using Range = FieldVector<double,dim>;
- // auto sourceTerm = [](const FieldVector<double,dim>& x){return Range{0.0, -1.0};};
-
Func2Tensor x3G_1 = [] (const Domain& x) {
return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
@@ -1816,7 +1926,6 @@ int main(int argc, char *argv[])
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
-
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
@@ -1824,7 +1933,6 @@ int main(int argc, char *argv[])
// log << basisContainer[0] << std::endl;
-
/////////////////////////////////////////////////////////
// Assemble the system
/////////////////////////////////////////////////////////
@@ -1867,8 +1975,7 @@ int main(int argc, char *argv[])
{
auto rowLocal = localView.tree().child(k).localIndex(arbitraryIndex);
row[k] = localView.index(rowLocal);
- std::cout << "row[k]:" << row[k] << std::endl;
-
+// std::cout << "row[k]:" << row[k] << std::endl;
}
///////////////////////////////////////////////////////////
@@ -1976,7 +2083,7 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////////////
- // Extract phi_alpha & M_alpha coefficients
+ // Extract phi_alpha & M_alpha coefficients
////////////////////////////////////////////////////////////////////////////////////
VectorCT phi_1, phi_2, phi_3;
@@ -2019,8 +2126,11 @@ int main(int argc, char *argv[])
printmatrix(std::cout, M_1, "Corrector-Matrix M_1", "--");
printmatrix(std::cout, M_2, "Corrector-Matrix M_2", "--");
printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
-
-
+ log << "Corrector-Matrix M_1: \n" << M_1 << std::endl;
+ log << " --------------------" << std::endl;
+ log << "Corrector-Matrix M_2: \n" << M_2 << std::endl;
+ log << " --------------------" << std::endl;
+ log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
@@ -2028,7 +2138,7 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////
// Substract Integral-mean
- //////////////////////////////////////////////////////////////////////////// // ERROR :
+ //////////////////////////////////////////////////////////////////////////// // ERROR
using SolutionRange = FieldVector<double,dim>;
std::cout << "check integralMean phi_1: " << std::endl;
@@ -2067,7 +2177,7 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
- // Write Solution in Logs
+ // Write Solution (Corrector Coefficients) in Logs
/////////////////////////////////////////////////////////
log << "\nSolution of Corrector problems:\n";
if(write_corrector_phi1)
@@ -2090,7 +2200,7 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
- //////////////////////////////////////////////////////////////////////////// // ERROR :
+ //////////////////////////////////////////////////////////////////////////// // ERROR
using SolutionRange = FieldVector<double,dim>;
auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
@@ -2160,7 +2270,7 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
- // Compute effective quantities
+ // Compute effective quantities: Elastic law & Prestrain
/////////////////////////////////////////////////////////
MatrixRT Q(0);
@@ -2190,7 +2300,8 @@ int main(int argc, char *argv[])
// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
}
printmatrix(std::cout, Q, "Matrix Q", "--");
-
+ log << "Computed Matrix Q: " << std::endl;
+ log << Q << std::endl;
// compute B_hat
for(size_t a = 0; a < 3; a++)
@@ -2203,7 +2314,8 @@ int main(int argc, char *argv[])
B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
}
-
+ log << "Computed prestrain B_hat: " << std::endl;
+ log << B_hat << std::endl;
@@ -2212,6 +2324,8 @@ int main(int argc, char *argv[])
Q.solve(Beff,B_hat);
+ log << "Computed prestrain B_eff: " << std::endl;
+ log << Beff << std::endl;
std::cout << "Beff : " << std::endl;
for(size_t i=0; i<3; i++)
@@ -2219,7 +2333,7 @@ int main(int argc, char *argv[])
std::cout <<"Beff[i]: " << Beff[i] << std::endl;
}
-
+ // TODO
// compute q1,q2,q3
FieldVector<double ,3> g1 = {1.0 , 0 , 0};
FieldVector<double ,3> g2 = {0 , 1.0 , 0};
@@ -2245,20 +2359,20 @@ int main(int argc, char *argv[])
std::cout << "Gamma: " << gamma << std::endl;
- log << "q1: " << q1_c << std::endl;
- log << "q2: " << q2_c << std::endl;
- log << "q3: " << q3_c << std::endl;
- log << "b1: " << Beff[1] << std::endl;
- log << "b2: " << Beff[2] << std::endl;
- log << "b3: " << Beff[3] << std::endl;
- log << "b1_hat: " << B_hat[1] << std::endl;
- log << "b2_hat: " << B_hat[2] << std::endl;
- log << "b3_hat: " << B_hat[3] << std::endl;
+ log << "computed q1: " << q1_c << std::endl;
+ log << "computed q2: " << q2_c << std::endl;
+ log << "computed q3: " << q3_c << std::endl;
+ log << "computed b1: " << Beff[0] << std::endl;
+ log << "computed b2: " << Beff[1] << std::endl;
+ log << "computed b3: " << Beff[2] << std::endl;
+ log << "computed b1_hat: " << B_hat[0] << std::endl;
+ log << "computed b2_hat: " << B_hat[1] << std::endl;
+ log << "computed b3_hat: " << B_hat[2] << std::endl;
-
+ // TODO
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
@@ -2279,6 +2393,13 @@ int main(int argc, char *argv[])
std::cout << "q1 : " << q1 << std::endl;
std::cout << "q2 : " << q2 << std::endl;
std::cout << "q3 should be between q1 and q2" << std::endl;
+ log << " --- analytic solutions: --- " << std::endl;
+ log << "b1 : " << b1 << std::endl;
+ log << "b2 : " << b2 << std::endl;
+ log << "b3 : " << b3 << std::endl;
+ log << "q1 : " << q1 << std::endl;
+ log << "q2 : " << q2 << std::endl;
+ log << "q3 should be between q1 and q2" << std::endl;
@@ -2291,7 +2412,7 @@ int main(int argc, char *argv[])
// };
auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
- return MatrixRT{{ (mu1*mu2/((theta*mu1 +(1-theta)*mu2)*muTerm(x)) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ return MatrixRT{{ ((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x)) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
auto zeroFunction = [] (const Domain& x) {
@@ -2300,7 +2421,7 @@ int main(int argc, char *argv[])
- FieldVector<double,3> testvector = {1.0/4.0 , 0.0 , 0.25};
+// FieldVector<double,3> testvector = {1.0/4.0 , 0.0 , 0.25};
// std::cout << "t[2]: " << testvector[2] << std::endl;
// std::cout << "muTerm(t):" << muTerm(testvector) << std::endl;
// auto Teest = symPhi_1_analytic(testvector);
@@ -2312,8 +2433,12 @@ int main(int argc, char *argv[])
auto L2error = computeL2Error(Basis_CE,phi_1,gamma,symPhi_1_analytic);
std::cout << "L2-Error: " << L2error << std::endl;
+ auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
+
+
auto L2Norm_Symphi = computeL2Error(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
- std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl; // TODO Not Needed
+ std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl; // TODO Not Needed
VectorCT zeroVec;
zeroVec.resize(Basis_CE.size());
@@ -2321,6 +2446,8 @@ int main(int argc, char *argv[])
auto L2Norm_analytic = computeL2Error(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
std::cout << "L2-Norm(analytic): " << L2Norm_analytic << std::endl;
+
+ log << "L2-Error (symmetric Gradient phi_1):" << L2errorTest << std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
@@ -2402,12 +2529,6 @@ int main(int argc, char *argv[])
VTKWriter<GridView> vtkWriter(gridView_CE);
-
-// vtkWriter.addVertexData(
-// correctorFunction_Test,
-// VTK::FieldInfo("corrector Test", VTK::FieldInfo::Type::vector, dim));
-
-
vtkWriter.addVertexData(
correctorFunction_1,
VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
@@ -2417,11 +2538,8 @@ int main(int argc, char *argv[])
vtkWriter.addVertexData(
correctorFunction_3,
VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
-
-
vtkWriter.write("CellProblem-result");
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
log.close();
-
}
--
GitLab