diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index f64fd0ba6735159275345fa67cdcc680cdcb2cac..0c1d7a605f8b24119440c4f89e626f91b887bc09 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -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