Skip to content
Snippets Groups Projects

WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis

Open WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis
feature/bendingIsometries into master
Open Klaus Böhnlein requested to merge feature/bendingIsometries into master
Viewing commit cfc28d5a
Prev Next
Show latest version
1 file
+ 145
− 97
Compare changes
  • Side-by-side
  • Inline
src/bending-isometries.cc
+ 145
− 97
@@ -54,13 +54,16 @@
// #include <dune/gfe/spaces/stiefelmatrix.hh>
#include <dune/gmsh4/gmsh4reader.hh>
#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
const int dim = 2;
using namespace Dune;
/**
* @brief Comparison of different function definitions for testing purposes.
* @brief Comparison of different function definitions for debugging purposes.
* TODO: remove later
*/
template<class LocalDiscreteKirchhoffFunction,class DeformationFunction, class GridView>
@@ -99,15 +102,14 @@ auto compare(LocalDiscreteKirchhoffFunction& kirchhoffFunction,
printvector(std::cout, defOut, "defOut:", "--");
printmatrix(std::cout, KirchhoffRot, "KirchhoffRot:", "--");
printmatrix(std::cout, defRot, "defRot:", "--");
printmatrix(std::cout, defRot*jacobian, "defRot*jacobian:", "--");
// printmatrix(std::cout, defRot*jacobian, "defRot*jacobian:", "--"); // deprecated: when ReducedCubicHermite returned values in globalCoordinates, this was actually neccessary to get the right values.
}
}
}
/**
* @brief Check isometry constraint at the nodes of the Triangulation
* @brief Test that checks the isometry constraint at the nodes of the Triangulation
*/
template<class LocalDiscreteKirchhoffFunction, class GridView>
auto nodewiseIsometryTest(LocalDiscreteKirchhoffFunction& deformationFunction,
@@ -148,15 +150,9 @@ auto nodewiseIsometryTest(DeformationFunction& deformationGridViewFunction)
// evaluate at nodes:
for (int i=0; i<geometry.corners(); i++)
{
// auto out = localDeformation(geometry.local(geometry.corner(i)));
/**
* @brief For the values of the derivative we need to multiply with the jacobian
* since ReducedCubicHermiteTriangle is not affine equivalent.
* compare with 'discreteglobalbasisfunction.hh' operator()
* in class 'DiscreteGlobalBasisFunctionDerivative'.
*/
auto rot = localRotation(geometry.local(geometry.corner(i)))*geometry.jacobian(geometry.corner(i));
auto rot = localRotation(geometry.local(geometry.corner(i)));
auto product = rot.transposed()*rot;
FieldMatrix<double,2,2> I = ScaledIdentityMatrix<double,2>(1);
auto diff = (rot.transposed()*rot) - I;
@@ -222,6 +218,10 @@ void interpolate(const Basis& basis, Coefficients& coefficients, Function& f)
// printvector(std::cout, value, "value:", "--");
// printmatrix(std::cout, der, "der:", "--");
// std::cout << "sin(x):" << sin(pos[0]) << std::endl;
// std::cout << "cos(x):" << cos(pos[0]) << std::endl;
// printvector(std::cout, pos, "pos: ", "--");
// Loop over components of power basis
for(std::size_t k=0; k<3 ; k++)
{
@@ -250,7 +250,7 @@ void interpolate(const Basis& basis, Coefficients& coefficients, Function& f)
/**
* @brief Data structure conversion from 'VectorSpaceCoefficients' to 'IsometryCoefficients'.
* TODO: Throw an error if rotation-part is not orthogonal?
* TODO: Throw an error if rotation-part is not orthogonal? (expensive to check)
*/
template<class DiscreteKirchhoffBasis, class CoefficientBasis, class Coefficients, class IsometryCoefficients>
void vectorToIsometryCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchhoffBasis,
@@ -294,7 +294,7 @@ void vectorToIsometryCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchh
auto globalIdx_2 = localView.index(localIdx_2);
// blockedInterleaved -- Version
currentDeformation[k] = vectorCoefficients[globalIdx_0[0]][globalIdx_0[1]] ;
currentDeformation[k] = vectorCoefficients[globalIdx_0[0]][globalIdx_0[1]];
currentDerivative_x[k] = vectorCoefficients[globalIdx_1[0]][globalIdx_1[1]];
currentDerivative_y[k] = vectorCoefficients[globalIdx_2[0]][globalIdx_2[1]];
@@ -308,13 +308,13 @@ void vectorToIsometryCoefficientMap(const DiscreteKirchhoffBasis& discreteKirchh
FieldMatrix<RT,3,3> rotMatrix(0);
for(std::size_t k=0; k<3 ; k++)
{
rotMatrix[0][k] = currentDerivative_x[k];
rotMatrix[1][k] = currentDerivative_y[k];
rotMatrix[2][k] = cross[k];
rotMatrix[k][0] = currentDerivative_x[k];
rotMatrix[k][1] = currentDerivative_y[k];
rotMatrix[k][2] = cross[k];
}
// printmatrix(std::cout, rotMatrix, "rotMatrix:", "--");
using namespace Dune::Indices;
isometryCoefficients[globalIdxOut][_1].set(rotMatrix); //TODO: Throw if matrix is not orthogonal?
isometryCoefficients[globalIdxOut][_1].set(rotMatrix); //TODO: Throw if matrix is not orthogonal?
isometryCoefficients[globalIdxOut][_0] = (RealTuple<RT, 3>)currentDeformation;
}
}
@@ -438,9 +438,14 @@ int main(int argc, char *argv[])
else
DUNE_THROW(Exception, "Unknown structured grid type '" << structuredGridType << "' found!");
} else {
std::cout << "Read GMSH grid." << std::endl;
std::string gridPath = parameterSet.get<std::string>("gridPath");
std::string gridFile = parameterSet.get<std::string>("gridFile");
grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(gridPath + "/" + gridFile));
GridFactory<GridType> factory;
Gmsh4::LagrangeGridCreator creator{factory};
Gmsh4Reader reader{creator};
reader.read(gridPath + "/" + gridFile);
grid = factory.createGrid();
}
const int numLevels = parameterSet.get<int>("numLevels");
@@ -462,10 +467,9 @@ int main(int argc, char *argv[])
using namespace Functions::BasisFactory;
auto deformationBasis = makeBasis(
gridView,
power<3>(reducedCubicHermiteTriangle(),
blockedInterleaved()));
auto deformationBasis = makeBasis(gridView,
power<3>(reducedCubicHermiteTriangle(),
blockedInterleaved()));
// The next basis is used to assign (nonlinear) degrees of freedom to the grid vertices.
@@ -474,11 +478,10 @@ int main(int argc, char *argv[])
CoefficientBasis coefficientBasis(gridView);
// A basis for the tangent space (used to set DirichletNodes)
auto tangentBasis = makeBasis(
gridView,
power<Coefficient::TangentVector::dimension>(
lagrange<1>(),
blockedInterleaved()));
auto tangentBasis = makeBasis(gridView,
power<Coefficient::TangentVector::dimension>(
lagrange<1>(),
blockedInterleaved()));
///////////////////////////////////////////
@@ -505,33 +508,28 @@ int main(int argc, char *argv[])
VectorSpaceCoefficients x(deformationBasis.size());
IdentityGridEmbedding<double> identityGridEmbedding; // TODO: Identity only for testing
IdentityGridEmbedding<double> identityGridEmbedding;
// interpolate(deformationBasis, x, identityGridEmbedding);
interpolate(deformationBasis, x, pythonInitialIterate); //TEST: interpolation on a python function
interpolate(deformationBasis, x, pythonInitialIterate); // Interpolation on a python function
// -------------------------------------------------------------------
///////////////////////////////////////////
// Test Stiefelmanifold-Class
//////////////////////////////////////////
std::cout << "Testing Stiefel-manifold class:" << std::endl;
// Stiefel<double,3,2> S;
// Dune::FieldMatrix<double,3,2> tangentVector = {{1.0, 0.0}, {0.0, 1.0}, {0.0, 0.0}};
// S.exp(S,tangentVector);
// std::cout << "Testing Stiefel-manifold class:" << std::endl;
// // Stiefel<double,3,2> S;
// // Dune::FieldMatrix<double,3,2> tangentVector = {{1.0, 0.0}, {0.0, 1.0}, {0.0, 0.0}};
// // S.exp(S,tangentVector);
Stiefel<double,3,2> S;
// Stiefel<double,3,2> S;
Dune::FieldMatrix<double,3,2> testMat(0);
// Dune::FieldMatrix<double,3,2> testMat(0);
S.matrix(testMat);
printmatrix(std::cout, testMat, "testMat", "--");
// S.matrix(testMat);
// printmatrix(std::cout, testMat, "testMat", "--");
std::cout << "Testing Stiefel-manifold class done." << std::endl;
// std::cout << "Testing Stiefel-manifold class done." << std::endl;
// exit(0);
// -------------------------------------------------------------------
@@ -553,23 +551,6 @@ int main(int argc, char *argv[])
IsometryCoefficients isometryCoefficients(coefficientBasis.size());
AIsometryCoefficients isometryCoefficients_adouble(coefficientBasis.size());
// Test: fill with dummy Values
// Dune::FieldVector<double, 7> dummyDeformationAndRotation({1,1,1,0,0,0,1});
// Coefficient dummyProduct(dummyDeformationAndRotation);
// for (auto &&v : vertices(gridView))
// isometryCoefficients[gridView.indexSet().index(v)] = dummyProduct;
// TEST: access product manifold components:
// Dune::index_constant<0> index_0;
// Dune::index_constant<1> index_1;
// auto deform = dummyProduct[index_0];
// auto rot = dummyProduct[index_1];
// FieldMatrix<double,3,3> rotMatrix(0);
// rot.matrix(rotMatrix);
// printvector(std::cout, deform.globalCoordinates(), "deform.globalCoordinates()", "--");
// printmatrix(std::cout, rotMatrix, "rotMatrix", "--");
std::cout << "Number of Elements in the grid:" << gridView.size(0)<< std::endl;
std::cout << "Number of Nodes in the grid:" << gridView.size(dim)<< std::endl;
std::cout << "deformationBasis.size():" << deformationBasis.size() << std::endl;
@@ -578,8 +559,11 @@ int main(int argc, char *argv[])
// isometryToVectorCoefficientMap(deformationBasis,coefficientBasis,x,isometryCoefficients);
// for(int i=0; i<x.size(); i++)
// std::cout<< "x[i]:" << x[i] << std::endl;
/**
* @brief Print coefficient vector:
*/
// for(int i=0; i<x.size(); i++)
// std::cout<< "x[i]:" << x[i] << std::endl;
/**
* @brief TEST: compute isometry error
@@ -587,21 +571,64 @@ int main(int argc, char *argv[])
auto initialDeformationFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationBasis, x);
// nodewiseIsometryTest(initialDeformationFunction);
//evaluate discrete function on the nodes:
// auto localInitialDeformationFunction = localFunction(initialDeformationFunction);
// auto localInitialRotationFunction = derivative(localInitialDeformationFunction);
// auto gridView = deformationGridViewFunction.basis().gridView();
// for (const auto &element : elements(gridView))
// {
// auto geometry = element.geometry();
// localInitialDeformationFunction.bind(element);
// localInitialRotationFunction.bind(element);
// // evaluate at nodes:
// for (int i=0; i<geometry.corners(); i++)
// {
// std::cout << "sin(x):" << sin(geometry.corner(i)[0]) << std::endl;
// std::cout << "cos(x):" << cos(geometry.corner(i)[0]) << std::endl;
// std::cout << "geometry.corner(i):" << geometry.corner(i) << std::endl;
// auto def = localInitialDeformationFunction(geometry.local(geometry.corner(i)));
// auto rot = localInitialRotationFunction(geometry.local(geometry.corner(i)))*geometry.jacobian(geometry.local(geometry.corner(i)));
// printvector(std::cout, def, "def: ", "--");
// printmatrix(std::cout, rot, "rot: ", "--");
// }
// }
// exit(0);
// Convert coefficient data structure from 'VectorSpaceCoefficients' to 'IsometryCoefficients'
vectorToIsometryCoefficientMap(deformationBasis,coefficientBasis,x,isometryCoefficients);
vectorToIsometryCoefficientMap(deformationBasis,coefficientBasis,x,isometryCoefficients_adouble);
// for(int i=0; i<isometryCoefficients.size(); i++)
// {
// using namespace Dune::Indices;
// // std::cout << "isometryCoefficients[i][_0].globalCoordinates(): " << isometryCoefficients[i][_0].globalCoordinates() << std::endl;
// printvector(std::cout, isometryCoefficients[i][_0].globalCoordinates(), "isometryCoefficients[i][_0].globalCoordinates(): ", "--");
// printvector(std::cout, isometryCoefficients[i][_1].director(0), "isometryCoefficients[i][_1].director(0): ", "--");
// printvector(std::cout, isometryCoefficients[i][_1].director(1), "isometryCoefficients[i][_1].director(1): ", "--");
// }
// exit(0);
// std::cout<< "isometryCoefficients[i].globalCoordinates():" << isometryCoefficients[i].globalCoordinates() << std::endl;
/**
* @brief TEST: conversion of coefficient vectors
* @brief TEST: conversion of coefficient vectors (back & forth)
*/
// VectorSpaceCoefficients x_convert = x;
// IsometryCoefficients isometryCoefficientsTMP(coefficientBasis.size());
// vectorToIsometryCoefficientMap(deformationBasis,coefficientBasis,x_convert,isometryCoefficientsTMP);
// isometryToVectorCoefficientMap(deformationBasis,coefficientBasis,x_convert,isometryCoefficientsTMP);
// //Check difference
// for(int i=0; i<x.size(); i++)
// std::cout<< "x[i]-x_convert[i]:" << x[i]-x_convert[i] << std::endl;
VectorSpaceCoefficients x_convert = x;
IsometryCoefficients isometryCoefficientsTMP(coefficientBasis.size());
vectorToIsometryCoefficientMap(deformationBasis,coefficientBasis,x_convert,isometryCoefficientsTMP);
isometryToVectorCoefficientMap(deformationBasis,coefficientBasis,x_convert,isometryCoefficientsTMP);
//Check difference
for(int i=0; i<x.size(); i++)
{
if((x[i]-x_convert[i]).two_norm() > 1e-10)
std::cout << "Coefficient conversion failed! with x[i]-x_convert[i]:" << x[i]-x_convert[i] << std::endl;
}
@@ -626,14 +653,46 @@ int main(int argc, char *argv[])
using DeformationBasis = decltype(deformationBasis);
// using LocalDKFunction = GFE::LocalDiscreteKirchhoffBendingIsometry<DeformationBasis, CoefficientBasis, IsometryCoefficients>;
using LocalDKFunction = GFE::LocalDiscreteKirchhoffBendingIsometry<DeformationBasis, CoefficientBasis, AIsometryCoefficients>;
LocalDKFunction localDKFunction(deformationBasis, coefficientBasis, isometryCoefficients_adouble);
// nodewiseIsometryTest(localDKFunction, gridView);
nodewiseIsometryTest(localDKFunction, gridView);
/**
* @brief Test compare different function definitions for the same coefficient vector:
* -function defined as powerBasis of reducedCubicHermite-Elements
* -localdiscretekirchhoffbendingisometry (constraint to be an isometry at the nodes)
*/
// compare(localDKFunction,initialDeformationFunction ,gridView);
// /**
// * @brief TEST : evaluate discrete function on the nodes:
// *
// */
// for (const auto &element : elements(gridView))
// {
// auto geometry = element.geometry();
// localDKFunction.bind(element);
// // Update local coefficients
// // localFunction_.updateLocalCoefficients(localSolution,element);
// // evaluate at nodes:
// for (int i=0; i<geometry.corners(); i++)
// {
// std::cout << "sin(x):" << sin(geometry.corner(i)[0]) << std::endl;
// std::cout << "cos(x):" << cos(geometry.corner(i)[0]) << std::endl;
// std::cout << "geometry.corner(i):" << geometry.corner(i) << std::endl;
// std::cout << "geometry.local(geometry.corner(i)):" << geometry.local(geometry.corner(i)) << std::endl;
// auto def = localDKFunction.evaluate(geometry.local(geometry.corner(i))).globalCoordinates();
// auto rot = localDKFunction.evaluateDerivative(geometry.local(geometry.corner(i)));
// printvector(std::cout, def, "def: ", "--");
// printmatrix(std::cout, rot, "rot: ", "--");
// }
// }
// exit(0);
// TEST: evaluation
// for (const auto &element : elements(gridView))
@@ -670,28 +729,27 @@ int main(int argc, char *argv[])
/**
* @brief Get Prestrain Tensor
* @brief Get effective prestrain Tensor
*/
Dune::FieldVector<adouble,3> Beffvec = parameterSet.get<Dune::FieldVector<adouble,3>>("effectivePrestrain", {0.0, 0.0, 0.0});
Dune::FieldMatrix<adouble,2,2> effectivePrestrain = {{(adouble)Beffvec[0], (adouble)Beffvec[2]}, {(adouble)Beffvec[2],(adouble)Beffvec[1]} };
printmatrix(std::cout, effectivePrestrain, "effective prestrain (Beff): ", "--");
// Dune::FieldMatrix<double,2,2> Beff(0);
// pyModule.get("effectivePrestrain").toC<Dune::FieldMatrix<double,2,2>>(Beff);
// parameterSet.get("effectivePrestrain").toC<FieldVector<double,3>>(Beffvec);
// Dune::FieldMatrix<double,2,2> Beff = { {Beffvec[0], Beffvec[1]}, {Beffvec[1], Beffvec[3]} };
// Assembler using ADOL-C
/**
* @brief Setup nonconforming energy
*/
auto localEnergy_nonconforming = std::make_shared<GFE::DiscreteKirchhoffBendingEnergy<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce);
// auto localEnergy_nonconforming = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyPrestrained<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce);
// auto localEnergy_conforming = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyConforming<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce);
//TEST
/**
* @brief Setup conforming energy (WIP)
*/
auto localEnergy_conforming = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyZienkiewicz<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce);
// auto localEnergy_conforming = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyZienkiewiczProjected<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce);
// Prestrain energy:
/**
* @brief Setup energy featuring prestrain
*/
auto localEnergy_prestrain = std::make_shared<GFE::DiscreteKirchhoffBendingEnergyPrestrained<CoefficientBasis, LocalDKFunction, decltype(localForce), ACoefficient>>(localDKFunction, localForce, effectivePrestrain);
LocalGeodesicFEADOLCStiffness<CoefficientBasis, Coefficient> localGFEADOLCStiffness_conforming(localEnergy_conforming.get());
@@ -726,14 +784,6 @@ int main(int argc, char *argv[])
isometryCoefficients,
dirichletNodes,
parameterSet);
// solver.setup(*grid,
// &assembler_nonconforming,
// isometryCoefficients,
// dirichletNodes,
// parameterSet.get<double>("tolerance",1e-12), /* tolerance */
// parameterSet.get<int>("maxProximalNewtonSteps",100), /* maxProximalNewtonSteps */
// parameterSet.get<double>("initialRegularization",1.0), /* initialRegularization */
// parameterSet.get<bool>("instrumented", 0)/* instrumented */);
// /////////////////////////////////////////////////////
// Solve!
@@ -748,7 +798,6 @@ int main(int argc, char *argv[])
////////////////////////////////
// Output result
////////////////////////////////
std::string baseNameDefault = "bending-isometries-";
std::string baseName = parameterSet.get("baseName", baseNameDefault);
std::string resultFileName = parameterSet.get("resultPath", "")
@@ -767,7 +816,6 @@ int main(int argc, char *argv[])
// VectorSpaceCoefficients identity(deformationBasis.dimension());
// auto foo = derivative(identityGridEmbedding);
//
// std::cout << "foo: \n" << foo({0,0}) << std::endl;
interpolate(deformationBasis, identity, identityGridEmbedding);
@@ -827,8 +875,8 @@ int main(int argc, char *argv[])
// }
// We need to subsample, because VTK cannot natively display real third-order functions
// SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
// SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, 3));
vtkWriter.write(resultFileName);
}