Skip to content
Snippets Groups Projects
Commit 7dc4dc22 authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

update

parent c7bc8b78
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/dune-microstructure.cc
+ 204
− 25
View file @ 7dc4dc22
......@@ -82,7 +82,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
// using Domain = typename LocalView::GridView::Codim<0>::Geometry::GlobalCoordinate;
// using FuncScalar = std::function< ScalarRT(const Domain&) >;
// using Func2Tensor = std::function< MatrixRT(const Domain&) >;
// using Func2Tensor = std::function< Matload_alpha1 ,rixRT(const Domain&) >;
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
......@@ -638,6 +638,9 @@ void assembleCellStiffness(const Basis& basis,
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
// std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
// std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
......@@ -716,7 +719,7 @@ void assembleCellLoad(const Basis& basis,
loadFunctional.bind(element);
const int localPhiOffset = localView.size();
std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
// BlockVector<double> elementRhs;
......@@ -739,6 +742,116 @@ void assembleCellLoad(const Basis& basis,
}
}
template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
auto energy(const Basis& basis,
LocalFunction1& muLocal,
LocalFunction2& lambdaLocal,
const MatrixFunction& matrixFieldFuncA,
const MatrixFunction& matrixFieldFuncB)
{
auto energy = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
auto localView = basis.localView();
auto matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
auto matrixFieldA = localFunction(matrixFieldAGVF);
auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
auto matrixFieldB = localFunction(matrixFieldBGVF);
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
for (const auto& e : elements(basis.gridView()))
{
localView.bind(e);
matrixFieldA.bind(e);
matrixFieldB.bind(e);
muLocal.bind(e);
lambdaLocal.bind(e);
FieldVector<double,1> elementEnergy(0);
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
int order = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
for (size_t pt=0; pt < quad.size(); pt++) {
const Domain& quadPos = quad[pt].position();
const double integrationElement = geometry.integrationElement(quadPos);
auto strain1 = matrixFieldA(quadPos);
// St.Venant-Kirchhoff stress
MatrixRT stressU(0);
stressU.axpy(2*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
double trace = 0;
for (int ii=0; ii<nCompo; ii++)
trace += strain1[ii][ii];
for (int ii=0; ii<nCompo; ii++)
stressU[ii][ii] += lambdaLocal(quadPos) * trace;
auto strain2 = matrixFieldB(quadPos);
// Local energy density: stress times strain
double energyDensity = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
energyDensity += stressU[ii][jj] * strain2[ii][jj];
elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
}
energy += elementEnergy;
}
return energy;
}
......@@ -1006,7 +1119,7 @@ auto integralMean(const Basis& basis,
area += quadPoint.weight() * integrationElement;
}
}
std::cout << "Domain-Area: " << area << std::endl;
// std::cout << "Domain-Area: " << area << std::endl;
for(const auto& element : elements(basis.gridView()))
......@@ -1070,7 +1183,7 @@ auto subtractIntegralMean(const Basis& basis,
for(size_t k=0; k<dim; k++)
{
std::cout << "Integral-Mean: " << IM[k] << std::endl;
// std::cout << "Integral-Mean: " << IM[k] << std::endl;
auto idx = childToIndexMap(basis,k);
for ( int i : idx)
......@@ -1290,6 +1403,30 @@ Func2Tensor x3G_3 = [] (const Domain& x) {
/////////////////////////////////////////////////////////////////////// TODO
// TODO : PrestrainImp.hh
double theta = 0.5;
double p1 = 1;
double p2 = 2;
using std::abs;
Func2Tensor B = [] (const Domain& x) {
double theta = 0.5;
double p1 = 1;
double p2 = 2;
if (abs(x[0])>= (theta/2) && x[2] > 0)
return MatrixRT{{p1, 0.0 , 0.0}, {0.0, p1, 0.0}, {0.0, 0.0, p1}};
if (abs(x[0]) < (theta/2) && x[2] < 0)
return MatrixRT{{p2, 0.0 , 0.0}, {0.0, p2, 0.0}, {0.0, 0.0, p2}};
else
return MatrixRT{{0.0, 0.0 , 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
//////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////
......@@ -1561,46 +1698,88 @@ printmatrix(std::cout, M_3, "Corrector-Matrix M_1", "--");
auto localPhi_1 = localFunction(correctorFunction_1);
auto localPhi_2 = localFunction(correctorFunction_2);
auto localPhi_3 = localFunction(correctorFunction_3);
subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
subtractIntegralMean(Basis_CE, localPhi_2 , phi_2);
subtractIntegralMean(Basis_CE, localPhi_3 , phi_3);
////////////////////////////////////////////////////////////////////////
// CHECK INTEGRAL-MEAN:
// auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
// Basis_CE,
// phi_1);
//
// auto AVPhi_1 = localFunction(AVFunction_1);
// auto A = integralMean(Basis_CE, AVPhi_1);
// for(size_t i=0; i<3; i++)
// {
// std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
// }
////////////////////////////////////////////////////
subtractIntegralMean(Basis_CE, localPhi_1 , phi_1);
// auto df = derivative(correctorFunction_1); // does not work :(
// INTEGRAL-MEAN TEST:
auto AVFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_1);
auto AVPhi_1 = localFunction(AVFunction_1);
auto A = integralMean(Basis_CE, AVPhi_1);
for(size_t i=0; i<3; i++)
const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
const std::array<decltype(localPhi_1)*, 3> phiContainer= {&localPhi_1, &localPhi_2, &localPhi_3};
const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
// compute Q
MatrixRT Q(0);
VectorCT tmp1,tmp2;
FieldVector<double,3> B_hat ;
for(size_t a = 0; a < 3; a++)
{
std::cout << "Integral-Mean (TEST) : " << A[i] << std::endl;
for(size_t b=0; b < 3; b++)
{
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, tmp1 ,x3MatrixBasis[b]); //
auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] );
Q[a][b] = coeffContainer[a]*tmp1 + GGterm; // seems symmetric... check positiv definitness?
}
}
printmatrix(std::cout, Q, "Matrix Q", "--");
for(size_t a = 0; a < 3; a++)
{
// FieldVector<double,3> ones = {1.0, 1.0, 1.0};
//
// phi_1 - IM*;
assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp2 ,B);
auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B );
// // const-function
// auto f1 = [](const auto& x)
// {
// return IM; // does not work..
// };
//
//
// VectorCT IM_1;
// interpolate(Basis_CE, IM_1, f1);
B_hat[a] = coeffContainer[a]*tmp2 + GGterm;
}
for(size_t i=0; i<3; i++)
{
std::cout << B_hat[i] << std::endl;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment