Skip to content
Snippets Groups Projects
Commit df015a54 authored by Lisa Julia Nebel's avatar Lisa Julia Nebel
Browse files

Clean up the usage of MIXED_SPACE in cosserat-continuum 

Now all combinations of dim = dimworld or dim != dimworld and MIXED_SPACE = 0 or MIXED_SPACE = 1 at least compile.
parent 0d09b494
No related branches found
No related tags found
1 merge request!72Cosserat-Continuum-Nonplanar
Hide whitespace changes
Inline Side-by-side
src/cosserat-continuum.cc
+ 243
− 188
View file @ df015a54
#define MIXED_SPACE 0
#include <config.h>
#include <fenv.h>
......@@ -37,7 +39,6 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
#include <dune/gfe/cosseratenergystiffness.hh>
......@@ -45,10 +46,16 @@
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/cosseratvtkreader.hh>
#include <dune/gfe/geodesicfeassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/mixedgfeassembler.hh>
#if MIXED_SPACE
#include <dune/gfe/mixedriemanniantrsolver.hh>
#else
#include <dune/gfe/geodesicfeassemblerwrapper.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/rigidbodymotion.hh>
#endif
#if HAVE_DUNE_VTK
#include <dune/vtk/vtkreader.hh>
......@@ -59,22 +66,17 @@
const int dim = 2;
const int dimworld = 2;
//#define MIXED_SPACE
// Order of the approximation space for the displacement
const int displacementOrder = 2;
// Order of the approximation space for the microrotations
const int rotationOrder = 2;
#ifndef MIXED_SPACE
#if !MIXED_SPACE
static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
// Image space of the geodesic fe functions
typedef RigidBodyMotion<double,3> TargetSpace;
// Tangent vector of the image space
const int blocksize = TargetSpace::TangentVector::dimension;
using TargetSpace = RigidBodyMotion<double, 3>;
#endif
using namespace Dune;
......@@ -84,6 +86,15 @@ int main (int argc, char *argv[]) try
// initialize MPI, finalize is done automatically on exit
Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
if (mpiHelper.rank()==0) {
std::cout << "DISPLACEMENTORDER = " << displacementOrder << ", ROTATIONORDER = " << rotationOrder << std::endl;
std::cout << "dim = " << dim << ", dimworld = " << dimworld << std::endl;
#if MIXED_SPACE
std::cout << "MIXED_SPACE = 1" << std::endl;
#else
std::cout << "MIXED_SPACE = 0" << std::endl;
#endif
}
// Start Python interpreter
Python::start();
Python::Reference main = Python::import("__main__");
......@@ -97,12 +108,8 @@ int main (int argc, char *argv[]) try
<< std::endl;
using namespace TypeTree::Indices;
#ifdef MIXED_SPACE
using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
std::vector<Rotation<double,3> > >;
#else
typedef std::vector<TargetSpace> SolutionType;
#endif
// parse data file
ParameterTree parameterSet;
......@@ -145,6 +152,8 @@ int main (int argc, char *argv[]) try
std::string structuredGridType = parameterSet["structuredGrid"];
if (structuredGridType == "cube") {
if (dim!=dimworld)
DUNE_THROW(GridError, "Please use FoamGrid and read in a grid for problems with dim != dimworld.");
lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
upper = parameterSet.get<FieldVector<double,dimworld> >("upper");
......@@ -183,7 +192,7 @@ int main (int argc, char *argv[]) try
GridView gridView = grid->leafGridView();
using namespace Dune::Functions::BasisFactory;
#ifdef MIXED_SPACE
const int dimRotation = Rotation<double,dim>::embeddedDim;
auto compositeBasis = makeBasis(
gridView,
......@@ -195,6 +204,13 @@ int main (int argc, char *argv[]) try
lagrange<rotationOrder>()
)
));
using CompositeBasis = decltype(compositeBasis);
auto deformationPowerBasis = makeBasis(
gridView,
power<3>(
lagrange<displacementOrder>()
));
typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
......@@ -202,10 +218,6 @@ int main (int argc, char *argv[]) try
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
#else
typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, displacementOrder> FEBasis;
FEBasis feBasis(gridView);
#endif
// /////////////////////////////////////////
// Read Dirichlet values
......@@ -261,43 +273,26 @@ int main (int argc, char *argv[]) try
if (orientationDirichletNodes[i][0])
for (int j=0; j<3; j++)
orientationDirichletDofs[i][j] = true;
#else
BitSetVector<1> dirichletNodes(feBasis.size(), false);
constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
BitSetVector<1> neumannNodes(feBasis.size(), false);
constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
for (size_t i=0; i<feBasis.size(); i++)
if (dirichletNodes[i][0])
for (int j=0; j<5; j++)
dirichletDofs[i][j] = true;
#endif
// //////////////////////////
// Initial iterate
// //////////////////////////
#ifdef MIXED_SPACE
SolutionType x;
x[_0].resize(deformationFEBasis.size());
x[_0].resize(compositeBasis.size({0}));
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
std::vector<FieldVector<double,3> > v;
Functions::interpolate(deformationFEBasis, v, pythonInitialDeformation);
Functions::interpolate(deformationPowerBasis, v, pythonInitialDeformation);
for (size_t i=0; i<x[_0].size(); i++)
x[_0][i] = v[i];
x[_1].resize(orientationFEBasis.size());
#else
SolutionType x(feBasis.size());
x[_1].resize(compositeBasis.size({1}));
#if !MIXED_SPACE
if (parameterSet.hasKey("startFromFile"))
{
std::shared_ptr<GridType> initialIterateGrid;
......@@ -314,35 +309,32 @@ int main (int argc, char *argv[]) try
initialIterateGrid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + initialIterateGridFilename));
}
SolutionType initialIterate;
std::vector<TargetSpace> initialIterate;
GFE::CosseratVTKReader::read(initialIterate, parameterSet.get<std::string>("initialIterateFilename"));
typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, 2> InitialBasis;
InitialBasis initialBasis(initialIterateGrid->leafGridView());
#ifdef PROJECTED_INTERPOLATION
using LocalInterpolationRule = LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
using LocalInterpolationRule = LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
#else
using LocalInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
using LocalInterpolationRule = LocalGeodesicFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
#endif
GFE::EmbeddedGlobalGFEFunction<InitialBasis,LocalInterpolationRule,TargetSpace> initialFunction(initialBasis,initialIterate);
auto powerBasis = makeBasis(
gridView,
power<7>(
lagrange<displacementOrder>()
));
std::vector<FieldVector<double,7> > v;
Dune::Functions::interpolate(feBasis,v,initialFunction);
DUNE_THROW(NotImplemented, "Replace scalar basis by power basis!");
for (size_t i=0; i<x.size(); i++)
x[i] = TargetSpace(v[i]);
//TODO: Interpolate does not work with an GFE:EmbeddedGlobalGFEFunction
//Dune::Functions::interpolate(powerBasis,v,initialFunction);
} else {
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
auto pythonInitialDeformation = Python::make_function<FieldVector<double,3> >(Python::evaluate(lambda));
std::vector<FieldVector<double,3> > v;
Functions::interpolate(feBasis, v, pythonInitialDeformation);
for (size_t i=0; i<x.size(); i++)
x[i].r = v[i];
for (size_t i=0; i<x.size(); i++) {
auto vTargetSpace = TargetSpace(v[i]);
x[_0][i] = vTargetSpace.r;
x[_1][i] = vTargetSpace.q;
}
}
#endif
......@@ -351,14 +343,30 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
#ifdef MIXED_SPACE
BlockVector<FieldVector<double,3> > identity(compositeBasis.size({0}));
Dune::Functions::interpolate(deformationPowerBasis, identity, [&](FieldVector<double,dimworld> x){ return x;});
BlockVector<FieldVector<double,3> > displacement(compositeBasis.size({0}));
if (dim == dimworld) {
CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
orientationFEBasis,x[_1],
resultPath + "mixed-cosserat_homotopy_0");
} else {
#if MIXED_SPACE
for (int i = 0; i < displacement.size(); i++) {
for (int j = 0; j < 3; j++)
displacement[i][j] = x[_0][i][j];
displacement[i] -= identity[i];
}
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationPowerBasis, displacement);
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, dimworld));
vtkWriter.write(resultPath + "mixed-cosserat_homotopy_0");
#else
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_0");
#endif
CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_0");
#endif
}
for (int i=0; i<numHomotopySteps; i++) {
double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
......@@ -396,83 +404,6 @@ int main (int argc, char *argv[]) try
materialParameters.report();
}
// Assembler using ADOL-C
#ifdef MIXED_SPACE
CosseratEnergyLocalStiffness<decltype(compositeBasis),
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
&neumannBoundary,
neumannFunction,
volumeLoad);
MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
RealTuple<double,3>,
Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
MixedGFEAssembler<decltype(compositeBasis),
RealTuple<double,3>,
Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
#else
std::shared_ptr<GFE::LocalEnergy<FEBasis,RigidBodyMotion<adouble,3> > > localCosseratEnergy;
if (dim==dimworld)
{
localCosseratEnergy = std::make_shared<CosseratEnergyLocalStiffness<FEBasis,3,adouble> >(materialParameters,
&neumannBoundary,
neumannFunction,
volumeLoad);
}
else
{
localCosseratEnergy = std::make_shared<NonplanarCosseratShellEnergy<FEBasis,3,adouble> >(materialParameters,
nullptr,
&neumannBoundary,
neumannFunction,
volumeLoad);
}
LocalGeodesicFEADOLCStiffness<FEBasis,
TargetSpace> localGFEADOLCStiffness(localCosseratEnergy.get());
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, localGFEADOLCStiffness);
#endif
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
#ifdef MIXED_SPACE
MixedRiemannianTrustRegionSolver<GridType,
decltype(compositeBasis),
DeformationFEBasis, RealTuple<double,3>,
OrientationFEBasis, Rotation<double,3> > solver;
#else
RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
#endif
solver.setup(*grid,
&assembler,
#ifdef MIXED_SPACE
deformationFEBasis,
orientationFEBasis,
#endif
x,
#ifdef MIXED_SPACE
deformationDirichletDofs,
orientationDirichletDofs,
#else
dirichletDofs,
#endif
tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
////////////////////////////////////////////////////////
// Set Dirichlet values
////////////////////////////////////////////////////////
......@@ -485,85 +416,208 @@ int main (int argc, char *argv[]) try
Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
// Extract object member functions as Dune functions
auto deformationDirichletValues = Python::make_function<FieldVector<double,3> >(dirichletValuesPythonObject.get("deformation"));
auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> >(dirichletValuesPythonObject.get("orientation"));
std::vector<FieldVector<double,3> > ddV;
std::vector<FieldMatrix<double,3,3> > dOV;
#ifdef MIXED_SPACE
Functions::interpolate(deformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
for (size_t j=0; j<x[_0].size(); j++)
if (deformationDirichletNodes[j][0])
x[_0][j] = ddV[j];
for (size_t j=0; j<x[_1].size(); j++)
if (orientationDirichletNodes[j][0])
x[_1][j].set(dOV[j]);
auto deformationDirichletValues = Python::make_function<FieldVector<double,3> > (dirichletValuesPythonObject.get("deformation"));
auto orientationDirichletValues = Python::make_function<FieldMatrix<double,3,3> > (dirichletValuesPythonObject.get("orientation"));
BlockVector<FieldVector<double,3> > ddV;
Dune::Functions::interpolate(deformationPowerBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
BlockVector<FieldMatrix<double,3,3> > dOV;
Dune::Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
for (int i = 0; i < compositeBasis.size({0}); i++) {
if (deformationDirichletDofs[i][0])
x[_0][i] = ddV[i];
}
for (int i = 0; i < compositeBasis.size({1}); i++)
if (orientationDirichletDofs[i][0])
x[_1][i].set(dOV[i]);
if (dim==dimworld) {
CosseratEnergyLocalStiffness<CompositeBasis, 3,adouble> localCosseratEnergy(materialParameters,
&neumannBoundary,
neumannFunction,
volumeLoad);
MixedLocalGFEADOLCStiffness<CompositeBasis,
RealTuple<double,3>,
Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy);
MixedGFEAssembler<CompositeBasis,
RealTuple<double,3>,
Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
#if MIXED_SPACE
MixedRiemannianTrustRegionSolver<GridType,
CompositeBasis,
DeformationFEBasis, RealTuple<double,3>,
OrientationFEBasis, Rotation<double,3> > solver;
solver.setup(*grid,
&mixedAssembler,
deformationFEBasis,
orientationFEBasis,
x,
deformationDirichletDofs,
orientationDirichletDofs, tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
#else
Functions::interpolate(feBasis, ddV, deformationDirichletValues, dirichletDofs);
Functions::interpolate(feBasis, dOV, orientationDirichletValues, dirichletDofs);
for (size_t j=0; j<x.size(); j++)
if (dirichletNodes[j][0])
{
x[j].r = ddV[j];
x[j].q.set(dOV[j]);
}
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
//Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
for (int i = 0; i < compositeBasis.size({0}); i++) {
for (int j = 0; j < 3; j ++) { // Displacement part
xTargetSpace[i].r[j] = x[_0][i][j];
dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
}
xTargetSpace[i].q = x[_1][i]; // Rotation part
for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
}
using GFEAssemblerWrapper = Dune::GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, TargetSpace, RealTuple<double, 3>, Rotation<double,3>>;
GFEAssemblerWrapper assembler(&mixedAssembler, deformationFEBasis);
RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace, GFEAssemblerWrapper> solver;
solver.setup(*grid,
&assembler,
xTargetSpace,
dirichletDofsTargetSpace,
tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
for (int i = 0; i < xTargetSpace.size(); i++) {
x[_0][i] = xTargetSpace[i].r;
x[_1][i] = xTargetSpace[i].q;
}
#endif
// /////////////////////////////////////////////////////
// Solve!
// /////////////////////////////////////////////////////
solver.setInitialIterate(x);
solver.solve();
x = solver.getSol();
} else {
#if MIXED_SPACE
DUNE_THROW(Exception, "Problems with dim != dimworld do not work for MIXED_SPACE = 1!");
#else
std::shared_ptr<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>> localGFEADOLCStiffness;
NonplanarCosseratShellEnergy<DeformationFEBasis, 3, adouble> localCosseratEnergyPlanar(materialParameters,
nullptr,
&neumannBoundary,
neumannFunction,
volumeLoad);
localGFEADOLCStiffness = std::make_shared<LocalGeodesicFEADOLCStiffness<DeformationFEBasis, TargetSpace>>(&localCosseratEnergyPlanar);
GeodesicFEAssembler<DeformationFEBasis,TargetSpace> assembler(gridView, *localGFEADOLCStiffness);
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
//Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
for (int i = 0; i < compositeBasis.size({0}); i++) {
for (int j = 0; j < 3; j ++) { // Displacement part
xTargetSpace[i].r[j] = x[_0][i][j];
dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
}
xTargetSpace[i].q = x[_1][i]; // Rotation part
for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
}
RiemannianTrustRegionSolver<DeformationFEBasis, TargetSpace> solver;
solver.setup(*grid,
&assembler,
xTargetSpace,
dirichletDofsTargetSpace,
tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
solver.setInitialIterate(xTargetSpace);
solver.solve();
xTargetSpace = solver.getSol();
for (int i = 0; i < xTargetSpace.size(); i++) {
x[_0][i] = xTargetSpace[i].r;
x[_1][i] = xTargetSpace[i].q;
}
#endif
}
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
#ifdef MIXED_SPACE
if (dim == dimworld) {
CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
orientationFEBasis,x[_1],
resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
#else
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_" + iAsAscii.str());
} else {
#if MIXED_SPACE
for (int i = 0; i < displacement.size(); i++) {
for (int j = 0; j < 3; j++) {
displacement[i][j] = x[_0][i][j];
}
displacement[i] -= identity[i];
}
auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3>>(deformationPowerBasis, displacement);
// We need to subsample, because VTK cannot natively display real second-order functions
SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(displacementOrder-1));
vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::scalar, 3));
vtkWriter.write(resultPath + "cosserat_homotopy_" + std::to_string(i+1));
#else
CosseratVTKWriter<GridType>::write<DeformationFEBasis>(deformationFEBasis, x, resultPath + "cosserat_homotopy_" + std::to_string(i+1));
#endif
}
}
// //////////////////////////////
// Output result
// //////////////////////////////
#ifndef MIXED_SPACE
#if !MIXED_SPACE
// Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
// This data may be used by other applications measuring the discretization error
BlockVector<TargetSpace::CoordinateType> xEmbedded(x.size());
for (size_t i=0; i<x.size(); i++)
xEmbedded[i] = x[i].globalCoordinates();
BlockVector<TargetSpace::CoordinateType> xEmbedded(compositeBasis.size({0}));
for (size_t i=0; i<compositeBasis.size({0}); i++)
for (size_t j=0; j<TargetSpace::TangentVector::dimension; j++)
xEmbedded[i][j] = j < 3 ? x[_0][i][j] : x[_1][i].globalCoordinates()[j-3];
std::ofstream outFile("cosserat-continuum-result-" + std::to_string(numLevels) + ".data", std::ios_base::binary);
MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
#endif
if (parameterSet.get<bool>("computeAverageNeumannDeflection", false) and parameterSet.hasKey("neumannValues")) {
// finally: compute the average deformation of the Neumann boundary
// That is what we need for the locking tests
FieldVector<double,3> averageDef(0);
#ifdef MIXED_SPACE
for (size_t i=0; i<x[_0].size(); i++)
if (neumannNodes[i][0])
averageDef += x[_0][i].globalCoordinates();
#else
for (size_t i=0; i<x.size(); i++)
if (neumannNodes[i][0])
averageDef += x[i].r;
#endif
averageDef /= neumannNodes.count();
if (mpiHelper.rank()==0 and parameterSet.hasKey("neumannValues"))
......@@ -571,6 +625,7 @@ int main (int argc, char *argv[]) try
std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << " "
<< ", average deflection: " << averageDef << std::endl;
}
}
} catch (Exception& e)
{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment