diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh index ab6f03ba85b2e1bb24bbaeb44839a28f9c45503e..4c8cba442cda47ac746cb1de0394013abc7178b4 100644 --- a/dune/microstructure/prestrainedMaterial.hh +++ b/dune/microstructure/prestrainedMaterial.hh @@ -11,6 +11,16 @@ #include <dune/microstructure/matrix_operations.hh> #include <dune/microstructure/voigthelper.hh> +//TEST dune-vtk +// #include <dune/vtk/function.hh> +#include <dune/vtk/vtkwriter.hh> +// #include <dune/vtk/vtkwriterbase.hh> +#include <dune/vtk/writers/unstructuredgridwriter.hh> +#include <dune/vtk/datacollectors/continuousdatacollector.hh> +#include <dune/vtk/datacollectors/discontinuousdatacollector.hh> +#include <dune/vtk/datacollectors/quadraticdatacollector.hh> +#include <dune/vtk/datacollectors/lagrangedatacollector.hh> + using namespace MatrixOperations; using namespace Dune::Functions; using namespace Dune::Functions::BasisFactory; @@ -563,7 +573,7 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase) std::cout << "wrote data to file:" + resultPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level) << std::endl; - //Version using the Dune::VTKSequenceWriter + //-- Version using the Dune::VTKSequenceWriter // std::shared_ptr<Dune::SubsamplingVTKWriter<GridView> > MaterialVtkWriter = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(gridView_, Dune::refinementLevels(subsamplingRefinement)); // MaterialVtkWriter->addCellData(indicatorFunction_, @@ -573,6 +583,17 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase) // Dune::VTKSequenceWriter<GridView> SequenceMaterialVTK(MaterialVtkWriter, resultPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level)); // SequenceMaterialVTK.write(vtkIndex); // std::cout << "wrote data to file:" + resultPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level) + "-" + std::to_string(vtkIndex) << std::endl; + + //--Version using dune-vtk + Dune::Vtk::DiscontinuousDataCollector<GridView> discontinuousDataCollector(gridView_); + // Dune::Vtk::YaspDataCollector<GridView, 4> yaspDataCollector(gridView_); + // Dune::Vtk::LagrangeDataCollector<GridView, 4> lagrangeDataCollector(gridView_); + Dune::Vtk::UnstructuredGridWriter writer(discontinuousDataCollector, Dune::Vtk::FormatTypes::ASCII, Dune::Vtk::DataTypes::FLOAT32); + + auto f2 = Dune::Functions::makeAnalyticGridViewFunction(indicatorFunction_ ,gridView_); + // writer.addPointData(f2 , "IndicatorFunction"); + writer.addCellData(f2 , "IndicatorFunction"); + writer.write(resultPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level) + "_DuneVTK"); return; }; diff --git a/src/macro-problem.cc b/src/macro-problem.cc index ef57b674b12dcdc305d34c59a52913ce335ac356..d184f3ebc00a3605581ffa997841e27a7a1625b7 100644 --- a/src/macro-problem.cc +++ b/src/macro-problem.cc @@ -25,6 +25,7 @@ #include <dune/grid/io/file/vtk.hh> #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh> +#include <dune/functions/gridfunctions/composedgridfunction.hh> #include <dune/functions/functionspacebases/lagrangebasis.hh> #include <dune/functions/functionspacebases/powerbasis.hh> #include <dune/functions/functionspacebases/interpolate.hh> @@ -65,8 +66,14 @@ //TEST dune-vtk -#include <dune/vtk/function.hh> +// #include <dune/vtk/function.hh> #include <dune/vtk/vtkwriter.hh> +// #include <dune/vtk/vtkwriterbase.hh> +#include <dune/vtk/writers/unstructuredgridwriter.hh> +#include <dune/vtk/datacollectors/continuousdatacollector.hh> +#include <dune/vtk/datacollectors/discontinuousdatacollector.hh> +#include <dune/vtk/datacollectors/quadraticdatacollector.hh> +#include <dune/vtk/datacollectors/lagrangedatacollector.hh> const int dim = 2; @@ -94,6 +101,46 @@ using namespace Dune; // #endif // } +struct Difference2 +{ + + double operator()(Dune::FieldMatrix<double, 3, 2> x, Dune::FieldMatrix<double, 2, 2> I) const + { + return ((x.transposed()*x)-I).frobenius_norm(); + } +}; + + +// Wrapper for global-coordinate functions F +template <class GridView, class F> +class GlobalKirchhoffFunction +{ + using Element = typename GridView::template Codim<0>::Entity; + using Geometry = typename Element::Geometry; + +public: + GlobalKirchhoffFunction (GridView const& gridView, F const& f) + : gridView_(gridView) + , f_(f) + {} + + void bind(Element const& element) { + geometry_.emplace(element.geometry()); + f_.bind(element); // we need to bind kirchhoff function to the element + } + void unbind() { geometry_.reset(); } + + auto operator() (typename Geometry::LocalCoordinate const& local) const + { + assert(!!geometry_); + return f_(geometry_->global(local)); + } + +private: + GridView gridView_; + F f_; + std::optional<Geometry> geometry_; +}; // Wrapper for global-coordinate functions F template <class GridView, class F> @@ -507,7 +554,20 @@ int main(int argc, char *argv[]) // //TEST dune-vtk: - VtkWriter writer{gridView, Dune::Vtk::FormatTypes::ASCII}; + + // Setup a DataCollector of order 4. + Dune::Vtk::LagrangeDataCollector<GridView, 4> lagrangeDataCollector(gridView); + Dune::Vtk::DiscontinuousDataCollector<GridView> discontinuousDataCollector(gridView); + + // VtkWriter writer{gridView, Dune::Vtk::FormatTypes::ASCII}; + // Dune::Vtk::VtkWriter writer{gridView, Dune::Vtk::FormatTypes::ASCII}; + + // Dune::Vtk::VtkWriterBase writer{gridView, Dune::Vtk::FormatTypes::ASCII}; + + + // Dune::Vtk::VtkWriterBase<GridView,Dune::Vtk::LagrangeDataCollector<GridView, 2> > writer{gridView, Dune::Vtk::FormatTypes::ASCII}; + // Vtk::UnstructuredGridWriter writer(lagrangeDataCollector, Vtk::FormatTypes::ASCII, Vtk::DataTypes::FLOAT32); + Vtk::UnstructuredGridWriter writer(discontinuousDataCollector, Vtk::FormatTypes::ASCII, Vtk::DataTypes::FLOAT32); /** @@ -520,21 +580,130 @@ int main(int argc, char *argv[]) flatLexicographic())); // auto normalBasis = makeBasis(gridView, // power<3>( - // lagrange<1>(), + // lagrang e<1>(), // blockedInterleaved())); - // TEST Compute + // TEST Compute isometry error function: + + + + +// // auto diff = (rot.transposed()*rot) - I; +// auto tmpFunction = GlobalKirchhoffFunction{gridView, [&](auto x) +// { +// return localDKFunctionDouble.evaluateDerivative(x); +// }}; + + + // auto isoErrorFunction = GlobalFunction{gridView, [&](auto x) + // { + // auto tmp = tmpFunction(x); + // auto out = (tmp.transposed()*tmp) - I; + // auto r = out.frobenius_norm(); + // return r; + + // }}; + + // auto isoErrorFunction = GlobalKirchhoffFunction{gridView, [&](auto x) + // // { + + // auto tmp = localDKFunctionDouble.evaluateDerivative(x); + + // auto out = (tmp.transposed()*tmp) - I; + + // auto r = out.frobenius_norm(); + + // return r; + + // }}; + + + + // writer.addPointData(isoErrorFunction, Dune::Vtk::FieldInfo{"isoErrorFunction dune-VTK",1, Vtk::RangeTypes::SCALAR}); + + + // ComposedGridFunction - approach: + auto deformationGradientFunction = derivative(deformationFunction); + auto deformationGradientLocalFunction = localFunction(deformationGradientFunction); + + + // auto normDiffFunction = [&](auto M) + // { + // return ((M.transposed()*M) - I).frobenius_norm(); + + // }; + auto normDiffFunction = [](Dune::FieldMatrix<double,3,2> M) -> double + { + FieldMatrix<double,2,2> Id = ScaledIdentityMatrix<double,2>(1); + return ((M.transposed()*M) - Id).frobenius_norm(); + }; + + // + + // auto gf_gridfunction = makeComposedGridFunction(normDiffFunction, std::ref(deformationGradientFunction)); + auto gf_gridfunction = makeComposedGridFunction(normDiffFunction, deformationGradientFunction); + // writer.addPointData(gf_gridfunction , Dune::Vtk::FieldInfo{"isoErrorFunction dune-VTK",1, Vtk::RangeTypes::SCALAR}); + // writer.addPointData(gf_gridfunction , VTK::FieldInfo("isoErrorFunction", VTK::FieldInfo::Type::scalar, 1)); + + vtkWriter.addVertexData(gf_gridfunction, VTK::FieldInfo("isoErrorFunction", VTK::FieldInfo::Type::scalar, 1)); + + + + auto identityF = Dune::Functions::makeAnalyticGridViewFunction([](auto&& x) -> Dune::FieldMatrix<double,2,2> { + return ScaledIdentityMatrix<double,2>(1); + + }, gridView); + // auto c = Dune::Function::makeComposedGridFunction(Difference2(),deformationGradientFunction,ScaledIdentityMatrix<double,2>(1)); + // auto c = Dune::Function::makeComposedGridFunction(Difference2(),deformationGradientFunction,ScaledIdentityMatrix<double,2>(1)); + + // auto isoErrorFunction =[](auto M, auto Id)-> double { + // return ((M.transposed()*M) - Id).frobenius_norm(); + // }; + + + // // auto c = makeComposedGridFunction(Difference2(),deformationGradientFunction,identityF); //works. + // auto c = makeComposedGridFunction(isoErrorFunction,deformationGradientFunction,identityF); + + FieldMatrix<double,2,2> I = ScaledIdentityMatrix<double,2>(1); + + auto isoErrorFunction =[&](auto M)-> double { + return ((M.transposed()*M) - I).frobenius_norm(); + }; + + + // auto c = makeComposedGridFunction(Difference2(),deformationGradientFunction,identityF); //works. + auto c = makeComposedGridFunction(isoErrorFunction,deformationGradientFunction); + + writer.addPointData(c, "IsometryError"); + // auto normalLambda = [deformationFunction](const FieldVector<double,2>& x) // { // // deformationfunction.derivative() ... //needs binding?! - - + // } + // auto isoErrorFunction = GlobalFunction{gridView, [&](auto x) + // { + // auto tmp = deformationGradientLocalFunction(x); + // auto out = (tmp.transposed()*tmp) - I; + // auto r = out.frobenius_norm(); + // return r; + // }}; + + // auto isoErrorFunction = GlobalFunction{gridView, [&](auto x) + // { + // return gf_gridfunction(x); + // }}; + + + // writer.addPointData(gf_gridfunction , VTK::FieldInfo("isoErrorFunction", VTK::FieldInfo::Type::scalar, 1)); + + // auto f1 = GlobalFunction{gridView, normDiff}; + // writer.addPointData(isoErrorFunction, VTK::FieldInfo("isoErrorFunction", VTK::FieldInfo::Type::scalar, 1)); /** @@ -557,8 +726,10 @@ int main(int argc, char *argv[]) // //TEST dune-vtk: // VtkWriter writer{gridView, Dune::Vtk::FormatTypes::ASCII}; - auto f2 = GlobalFunction{gridView, [&](auto x) {return pythonAnalyticalSolution(x); }}; - writer.addPointData(f2, Dune::Vtk::FieldInfo{"displacement_analytical dune-VTK",3, Vtk::RangeTypes::VECTOR}); + // auto f2 = GlobalFunction{gridView, [&](auto x) {return pythonAnalyticalSolution(x); }}; + auto f2 = Dune::Functions::makeAnalyticGridViewFunction( pythonAnalyticalSolution ,gridView); + // writer.addPointData(f2, Dune::Vtk::FieldInfo{"displacement_analytical dune-VTK",3, Vtk::RangeTypes::VECTOR}); + writer.addPointData(f2, "pythonAnalytic"); @@ -595,8 +766,10 @@ int main(int argc, char *argv[]) //dune-vtk - auto f3 = GlobalFunction{gridView, [&](auto x) {return pythonSurfaceNormal(x); }}; - writer.addPointData(f3, Dune::Vtk::FieldInfo{"surfaceNormal_analytical dune-VTK",3, Vtk::RangeTypes::VECTOR}); + // auto f3 = GlobalFunction{gridView, [&](auto x) {return pythonSurfaceNormal(x); }}; + auto f3 = Dune::Functions::makeAnalyticGridViewFunction(pythonSurfaceNormal ,gridView); + // writer.addPointData(f3, Dune::Vtk::FieldInfo{"surfaceNormal_analytical dune-VTK",3, Vtk::RangeTypes::VECTOR}); + writer.addPointData(f3, "pythonSurfaceNormal"); } } @@ -616,8 +789,10 @@ int main(int argc, char *argv[]) - writer.addPointData(displacementFunction, Dune::Vtk::FieldInfo{"displacement dune-VTK",3, Vtk::RangeTypes::VECTOR}); - writer.addPointData(surfaceNormalDiscrete, Dune::Vtk::FieldInfo{"surfaceNormal_discrete dune-VTK",3, Vtk::RangeTypes::VECTOR}); + // writer.addPointData(displacementFunction, Dune::Vtk::FieldInfo{"displacement dune-VTK",3, Vtk::RangeTypes::VECTOR}); + // writer.addPointData(surfaceNormalDiscrete, Dune::Vtk::FieldInfo{"surfaceNormal_discrete dune-VTK",3, Vtk::RangeTypes::VECTOR}); + writer.addPointData(displacementFunction, "displacement dune-VTK"); + writer.addPointData(surfaceNormalDiscrete, "surfaceNormal_discrete dune-VTK"); // dune-vtk