Skip to content
Snippets Groups Projects
Commit 31cb5912 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Fix writing of periodic deformations

parent f1a700a7
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/quasiconvexity-test.cc
+ 102
− 110
View file @ 31cb5912
......@@ -185,7 +185,107 @@ struct VossDensity final
double c_;
};
template <class Basis, class Vector>
void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldMatrix<double,2,2>& F0, std::string filename)
{
// Construct basis
using GridView = typename Basis::GridView;
using NonperiodicBasis = Functions::LagrangeBasis<GridView, order>;
NonperiodicBasis nonperiodicBasis(basis.gridView());
using namespace Dune::Functions::BasisFactory;
auto nonperiodicPowerBasis = makeBasis(
basis.gridView(),
power<dim>(
lagrange<order>(),
blockedInterleaved()
));
// Add the underlying affine displacement. For this we need
// a representation in a non-periodic basis.
BCRSMatrix<double> matrix;
assembleGlobalBasisTransferMatrix(matrix, basis, nonperiodicBasis);
Vector nonperiodicX(nonperiodicBasis.size());
matrix.mv(periodicX, nonperiodicX);
Vector affineDisplacement;
Dune::Functions::interpolate(nonperiodicPowerBasis, affineDisplacement, [&F0](FieldVector<double,dim> pos)
{
return FieldVector<double,2>{ (F0[0][0]-1)*pos[0] + F0[0][1]*pos[1], F0[1][0]*pos[0] + (F0[1][1]-1)*pos[1] };
});
Vector displacement = nonperiodicX;
displacement += affineDisplacement;
// Compute the determinant per element, for better understanding of the SolutionType
Functions::LagrangeBasis<GridView,0> p0Basis(basis.gridView());
std::vector<double> deformationDeterminant(p0Basis.size());
// K = 0.5 F.frobenius_norm2() / F.determinant();
std::vector<double> K(p0Basis.size());
auto localView = nonperiodicBasis.localView();
auto localP0View = p0Basis.localView();
for (const auto& element : elements(basis.gridView()))
{
localView.bind(element);
localP0View.bind(element);
const auto& localFiniteElement = localView.tree().finiteElement();
Vector localConfiguration(localView.size());
for (size_t i=0; i<localConfiguration.size(); i++)
localConfiguration[i] = nonperiodicX[localView.index(i)];
// store gradients of shape functions and base functions
std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients(localFiniteElement.size());
std::vector<Dune::FieldVector<double,dim> > gradients(localFiniteElement.size());
auto p = element.geometry().center();
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(p);
// get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(p, referenceGradients);
// compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
Dune::FieldMatrix<double,dim,dim> derivative(0);
for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<dim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
derivative += F0;
//std::cout << derivative << std::endl;
deformationDeterminant[localP0View.index(0)] = derivative.determinant();
K[localP0View.index(0)] = 0.5 * derivative.frobenius_norm2() / derivative.determinant();
}
std::cout << "K range: " << (*std::min_element(K.begin(), K.end()))
<< ", " << (*std::max_element(K.begin(), K.end())) << std::endl;
auto deformationDeterminantFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, deformationDeterminant);
auto KFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, K);
// Actually write the displacement function
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(nonperiodicBasis, displacement);
//auto localDisplacementFunction = localFunction(displacementFunction);
VTKWriter<GridView> vtkWriter(basis.gridView());
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
//vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addCellData(deformationDeterminantFunction,
VTK::FieldInfo("determinant", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.addCellData(KFunction,
VTK::FieldInfo("K", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write(filename);
}
int main (int argc, char *argv[]) try
{
......@@ -486,22 +586,7 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
SolutionType identity;
Dune::Functions::interpolate(powerBasis, identity, [&F0](FieldVector<double,dim> pos)
{
return FieldVector<double,2>{ (F0[0][0]-1)*pos[0] + F0[0][1]*pos[1], F0[1][0]*pos[0] + (F0[1][1]-1)*pos[1] };
});
SolutionType displacement = x;
displacement += identity;
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
auto localDisplacementFunction = localFunction(displacementFunction);
// Write the initial iterate
VTKWriter<GridView> vtkWriter(gridView);
vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(resultPath + "quasiconvexity-test_homotopy_0");
writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test_homotopy_0");
for (int i=0; i<numHomotopySteps; i++)
{
......@@ -695,101 +780,8 @@ int main (int argc, char *argv[]) try
std::cout << "Energy part 2: " << assembler.computeEnergy(x) << std::endl;
#endif
// Compute the displacement
auto displacement = x;
displacement += identity;
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
// Compute the determinant per element, for better understanding of the SolutionType
Dune::Functions::LagrangeBasis<GridView,0> p0Basis(gridView);
std::vector<double> deformationDeterminant(p0Basis.size());
// K = 0.5 F.frobenius_norm2() / F.determinant();
std::vector<double> K(p0Basis.size());
#if 0
auto deformationFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, x);
#if 1
auto deformationGradient = derivative(deformationFunction);
auto localDeformationGradient = localFunction(deformationGradient);
for (const auto& element : elements(gridView))
{
localDeformationGradient.bind(element);
auto F = localDeformationGradient(element.geometry().center());
std::cout << F << std::endl;
}
#else
auto localDeformation = localFunction(deformationFunction);
for (const auto& element : elements(gridView))
{
localDeformation.bind(element);
auto localDeformationGradient = derivative(localDeformation);
auto F = localDeformationGradient(element.geometry().center());
std::cout << F << std::endl;
}
#endif
#else
auto localView = feBasis.localView();
auto localP0View = p0Basis.localView();
for (const auto& element : elements(gridView))
{
localView.bind(element);
localP0View.bind(element);
const auto& localFiniteElement = localView.tree().finiteElement();
SolutionType localConfiguration(localView.size());
for (size_t i=0; i<localConfiguration.size(); i++)
localConfiguration[i] = x[localView.index(i)];
// store gradients of shape functions and base functions
std::vector<Dune::FieldMatrix<double,1,dim> > referenceGradients(localFiniteElement.size());
std::vector<Dune::FieldVector<double,dim> > gradients(localFiniteElement.size());
auto p = element.geometry().center();
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(p);
// get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(p, referenceGradients);
// compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
Dune::FieldMatrix<double,dim,dim> derivative(0);
for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<dim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
derivative += F0;
//std::cout << derivative << std::endl;
deformationDeterminant[localP0View.index(0)] = derivative.determinant();
K[localP0View.index(0)] = 0.5 * derivative.frobenius_norm2() / derivative.determinant();
}
#endif
std::cout << "p: " << parameterSet.sub("materialParameters")["p"] << ", "
<< "K range: " << (*std::min_element(K.begin(), K.end()))
<< ", " << (*std::max_element(K.begin(), K.end())) << std::endl;
auto deformationDeterminantFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, deformationDeterminant);
auto KFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, K);
VTKWriter<GridView> vtkWriter(gridView);
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
vtkWriter.addCellData(deformationDeterminantFunction,
VTK::FieldInfo("determinant", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.addCellData(KFunction,
VTK::FieldInfo("K", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write(resultPath + "quasiconvexity-test_homotopy_" + std::to_string(i+1));
writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test_homotopy_" + std::to_string(i+1));
}
// Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment