diff --git a/test/localgeodesicfestiffnesstest.cc b/test/localgeodesicfestiffnesstest.cc index 4724473fac3d9721328bc4854d3e76d7756eb189..dbaf6ee68fbde875e0351cdc9ae295d945c333eb 100644 --- a/test/localgeodesicfestiffnesstest.cc +++ b/test/localgeodesicfestiffnesstest.cc @@ -21,7 +21,7 @@ using namespace Dune; - +/** \brief A special energy functional of which I happen to be able to compute the Hessian */ template<class GridView, class TargetSpace> class TestEnergyLocalStiffness : public LocalGeodesicFEStiffness<GridView,TargetSpace> @@ -43,15 +43,6 @@ public: RT energy (const Entity& e, const std::vector<TargetSpace>& localSolution) const; - TestEnergyLocalStiffness(TargetSpace anchor, int idx) - : anchor_(anchor), - idx_(idx) - {} - - TargetSpace anchor_; - - int idx_; - }; template <class GridView, class TargetSpace> @@ -59,7 +50,8 @@ typename TestEnergyLocalStiffness<GridView, TargetSpace>::RT TestEnergyLocalStif energy(const Entity& element, const std::vector<TargetSpace>& localSolution) const { - return TargetSpace::distance(anchor_, localSolution[idx_]) * TargetSpace::distance(anchor_, localSolution[idx_]); + return TargetSpace::distance(localSolution[0], localSolution[1]) + * TargetSpace::distance(localSolution[0], localSolution[1]); } @@ -111,10 +103,7 @@ void testUnitVector3d() int nTestPoints = testPoints.size(); - TargetSpace anchor = testPoints[0]; - int idx = 0; - - TestEnergyLocalStiffness<typename GridType::LeafGridView, TargetSpace> assembler(anchor, idx); + TestEnergyLocalStiffness<typename GridType::LeafGridView, TargetSpace> assembler; // Set up elements of S^2 std::vector<TargetSpace> coefficients(domainDim+1); @@ -136,8 +125,13 @@ void testUnitVector3d() std::cout << "fdHessian:\n"; printmatrix(std::cout, fdHessian, "fdHessian", "--"); - FieldMatrix<double,3,3> embeddedHessian = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(anchor, coefficients[idx]); - FieldMatrix<double,2,2> hessian; + Matrix<FieldMatrix<double,3,3> > embeddedHessian(nDofs,nDofs); + + for (size_t j=0; j<nDofs; j++) + for (size_t k=0; k<nDofs; k++) + embeddedHessian[j][k] = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients[j], + coefficients[j]); + Matrix<FieldMatrix<double,2,2> > hessian(nDofs,nDofs); // transform to local tangent space bases const int blocksize = 2; @@ -154,7 +148,6 @@ void testUnitVector3d() } -#if 0 for (size_t j=0; j<nDofs; j++) for (size_t k=0; k<nDofs; k++) { @@ -163,14 +156,9 @@ void testUnitVector3d() hessian[j][k] = tmp.rightmultiplyany(orthonormalFramesTransposed[k]); } -#else - Dune::FieldMatrix<double,blocksize,embeddedBlocksize> tmp; - Dune::FMatrixHelp::multMatrix(orthonormalFrames[idx],embeddedHessian,tmp); - hessian = tmp.rightmultiplyany(orthonormalFramesTransposed[idx]); - -#endif - std::cout << "hessian:\n" << hessian << std::endl; + std::cout << "hessian:" << std::endl; + printmatrix(std::cout, hessian, "hessian", "--"); }