From 10d72af37157c511b16e5d68a9e85d0ac9d7a42b Mon Sep 17 00:00:00 2001
From: "felix.albrecht" <felix.albrecht@4028485c-44d9-4cde-a312-5a4635ee2db9>
Date: Wed, 6 Apr 2011 18:00:46 +0000
Subject: [PATCH] added example, fixed Makefiles and CmakeLists.txt
git-svn-id: https://dune.mathematik.uni-freiburg.de/svn/dune-fem-functionals/trunk@114 4028485c-44d9-4cde-a312-5a4635ee2db9
---
.gitignore | 9 +
CMakeLists.txt | 4 +
configure.ac | 1 +
dune/fem/functional/Makefile.am | 4 +-
examples/Makefile.am | 2 +-
examples/finite_element/main.cc | 6 -
.../Makefile.am | 35 ++
.../finite_element_local_operation.cc | 315 ++++++++++++++++++
.../macrogrids/unitcube1.dgf | 18 +
.../macrogrids/unitcube2.dgf | 18 +
.../macrogrids/unitcube3.dgf | 18 +
11 files changed, 421 insertions(+), 9 deletions(-)
create mode 100644 examples/finite_element_local_operation/Makefile.am
create mode 100644 examples/finite_element_local_operation/finite_element_local_operation.cc
create mode 100644 examples/finite_element_local_operation/macrogrids/unitcube1.dgf
create mode 100644 examples/finite_element_local_operation/macrogrids/unitcube2.dgf
create mode 100644 examples/finite_element_local_operation/macrogrids/unitcube3.dgf
diff --git a/.gitignore b/.gitignore
index ccf1d3dab..5b3ff7c36 100644
--- a/.gitignore
+++ b/.gitignore
@@ -62,3 +62,12 @@ examples/finite_element/test
examples/finite_element/test-main.o
tests/common/.deps/
tests/operator/.deps/
+CMakeLists.txt.user.2.1pre1
+CMakeLists.txt.user.2.3pre1
+examples/finite_element_local_operation/.deps
+examples/finite_element_local_operation/finite_element_local_operation
+examples/finite_element_local_operation/finite_element_local_operation-finite_element_local_operation.o
+tests/common/localbasefunction_test
+tests/common/localbasefunction_test-localbasefunction_test.o
+tests/operator/ellipticfiniteelementoperator_test
+tests/operator/ellipticfiniteelementoperator_test-ellipticfiniteelementoperator_test.o
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c311d8a3f..53f5a3e1e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -129,4 +129,8 @@ TARGET_LINK_LIBRARIES( localbasefunction_test ${COMMON_LIBS} )
ADD_EXECUTABLE( ellipticfiniteelement_example "examples/finite_element/main.cc" ${common} ${grid} ${istl} ${fem} ${femhowto} ${femtools} ${femfunctionals} )
TARGET_LINK_LIBRARIES( ellipticfiniteelement_example ${COMMON_LIBS} )
+ADD_EXECUTABLE( ellipticfiniteelement_local_operation_example "examples/finite_element_local_operation/finite_element_local_operation.cc" ${common} ${grid} ${istl} ${fem} ${femhowto} ${femtools} ${femfunctionals} )
+TARGET_LINK_LIBRARIES( ellipticfiniteelement_local_operation_example ${COMMON_LIBS} )
+
+
#HEADERCHECK( ${header} )
diff --git a/configure.ac b/configure.ac
index 5b502b761..2910eed3b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -35,6 +35,7 @@ AC_CONFIG_FILES([
examples/Makefile
examples/poisson_new/Makefile
examples/finite_element/Makefile
+ examples/finite_element_local_operation/Makefile
tests/Makefile
tests/common/Makefile
tests/operator/Makefile
diff --git a/dune/fem/functional/Makefile.am b/dune/fem/functional/Makefile.am
index 5634c18af..4288688d4 100644
--- a/dune/fem/functional/Makefile.am
+++ b/dune/fem/functional/Makefile.am
@@ -1,4 +1,4 @@
-ltwofunctionalincludedir = $(includedir)/dune/fem/functional
-ltwofunctionalinclude_HEADERS = ltwo.hh
+finiteelementincludedir = $(includedir)/dune/fem/functional
+finiteelementinclude_HEADERS = finiteelement.hh
include $(top_srcdir)/am/global-rules
diff --git a/examples/Makefile.am b/examples/Makefile.am
index f93362081..72566d903 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -1,6 +1,6 @@
# $Id$
-SUBDIRS = poisson_new finite_element
+SUBDIRS = poisson_new finite_element finite_element_local_operation
if BUILD_DOCS
# TODO: set up documentation tree automatically
diff --git a/examples/finite_element/main.cc b/examples/finite_element/main.cc
index c217161b2..9092f8a9c 100644
--- a/examples/finite_element/main.cc
+++ b/examples/finite_element/main.cc
@@ -144,24 +144,18 @@ public:
{
ret = 1.0;
- std::cout << "evaluate at ";
-
if (!(arg[0] > 0.0)) {
// bottom
- std::cout << "bottom" << std::endl;
ret = 1.0;
} else if (!(arg[0] < 1.0)) {
// top
- std::cout << "top" << std::endl;
ret = 1.0;
}
if (!(arg[1] > 0.0)) {
// left
- std::cout << "left" << std::endl;
ret = 1.0;
} else if (!(arg[1] < 1.0)) {
// right
- std::cout << "right" << std::endl;
ret = 1.0;
}
}
diff --git a/examples/finite_element_local_operation/Makefile.am b/examples/finite_element_local_operation/Makefile.am
new file mode 100644
index 000000000..9b727f1f2
--- /dev/null
+++ b/examples/finite_element_local_operation/Makefile.am
@@ -0,0 +1,35 @@
+GRIDDIM=2
+GRIDTYPE=YASPGRID
+POLORDER=1
+
+noinst_PROGRAMS = finite_element_local_operation
+
+finite_element_local_operation_SOURCES = finite_element_local_operation.cc
+
+finite_element_local_operation_CPPFLAGS = $(AM_CPPFLAGS) \
+ $(DUNEMPICPPFLAGS) \
+ $(UG_CPPFLAGS) \
+ $(AMIRAMESH_CPPFLAGS) \
+ $(ALBERTA_CPPFLAGS) \
+ $(ALUGRID_CPPFLAGS) -DGRIDDIM=$(GRIDDIM) -D$(GRIDTYPE) -DPOLORDER=$(POLORDER) -Wall -O0 -DDEBUG -funroll-loops -g -ggdb -fno-strict-aliasing -std=c++0x
+# The libraries have to be given in reverse order (most basic libraries
+# last). Also, due to some misunderstanding, a lot of libraries include the
+# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
+# here as well.
+finite_element_local_operation_LDADD = \
+ $(DUNE_LDFLAGS) $(DUNE_LIBS) \
+ $(ALUGRID_LDFLAGS) $(ALUGRID_LIBS) \
+ $(ALBERTA_LDFLAGS) $(ALBERTA_LIBS) \
+ $(AMIRAMESH_LDFLAGS) $(AMIRAMESH_LIBS) \
+ $(UG_LDFLAGS) $(UG_LIBS) \
+ $(DUNEMPILIBS) \
+ $(LDADD)
+finite_element_local_operation_LDFLAGS = $(AM_LDFLAGS) \
+ $(DUNEMPILDFLAGS) \
+ $(UG_LDFLAGS) \
+ $(AMIRAMESH_LDFLAGS) \
+ $(ALBERTA_LDFLAGS) \
+ $(ALUGRID_LDFLAGS) \
+ $(DUNE_LDFLAGS)
+
+include $(top_srcdir)/am/global-rules
diff --git a/examples/finite_element_local_operation/finite_element_local_operation.cc b/examples/finite_element_local_operation/finite_element_local_operation.cc
new file mode 100644
index 000000000..9092f8a9c
--- /dev/null
+++ b/examples/finite_element_local_operation/finite_element_local_operation.cc
@@ -0,0 +1,315 @@
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <iostream>
+#include <vector>
+
+// disable warnings about problems in dune headers
+#include <dune/fem-tools/header/disablewarnings.hh>
+
+// dune-common includes
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+
+// dune-grid includes
+#include <dune/grid/utility/gridtype.hh>
+
+// dune-istl includes
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/bvector.hh>
+
+// dune-fem includes
+#include <dune/fem/misc/mpimanager.hh>
+#include <dune/fem/gridpart/gridpart.hh>
+#include <dune/fem/space/lagrangespace.hh>
+#include <dune/fem/storage/vector.hh>
+#include <dune/fem/function/adaptivefunction/adaptivefunction.hh>
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+
+// reenable warnings about problems in dune headers
+#include <dune/fem-tools/header/enablewarnings.hh>
+
+// dune-fem-functionals includes
+#include <dune/fem/constraints/dirichlet.hh>
+#include <dune/fem/subspace/subspaces.hh>
+#include <dune/fem/operator/finiteelement.hh>
+#include <dune/fem/functional/finiteelement.hh>
+#include <dune/fem/container/factory.hh>
+#include <dune/fem/solver/femassembler.hh>
+
+// dune-fem-tools includes
+#include <dune/fem-tools/function/functiontools.hh>
+#include <dune/fem-tools/space/projection.hh>
+
+using namespace Dune::Functionals;
+
+#ifndef POLORDER
+const int polOrder = 1;
+#else
+const int polOrder = POLORDER;
+#endif
+
+/**
+ * \brief Represents the elliptic operation a(x) \gradient u(x) \gradient v(x) for given u, v, x.
+ * In this case, a = 1.
+ **/
+class EllipticOperation
+{
+public:
+ template <class FirstLocalFunctionType, class SecondLocalFunctionType, class LocalPointType>
+ double operate(const FirstLocalFunctionType& firstLocalFunction, const SecondLocalFunctionType& secondLocalFunction,
+ const LocalPointType& localPoint) const
+ {
+ // init return value
+ double ret = 0.0;
+
+ // some types we will need
+ typedef typename SecondLocalFunctionType::EntityType EntityType;
+
+ typedef typename FirstLocalFunctionType::JacobianRangeType JacobianRangeType;
+
+ // first gradient
+ JacobianRangeType firstGradient(0.0);
+ firstLocalFunction.jacobian(localPoint, firstGradient);
+
+ // second gradient
+ JacobianRangeType secondGradient(0.0);
+ secondLocalFunction.jacobian(localPoint, secondGradient);
+
+ const double product = firstGradient[0] * secondGradient[0];
+
+ // 1.0 * \gradient u(x) \gradient v(x)
+ ret = 1.0 * product;
+
+ // return
+ return ret;
+ }
+
+}; // end class EllipticOperation
+
+/**
+ * \brief Represents the product operation f(x) v(x) for given v, x.
+ * In this case, f = 1.
+ **/
+class ProductOperation
+{
+public:
+ template <class LocalFunctionType, class LocalPointType>
+ double operate(const LocalFunctionType& localFunction, const LocalPointType& localPoint) const
+ {
+ // init return value
+ double ret = 0.0;
+
+ // some types we will need
+ typedef typename LocalFunctionType::RangeType RangeType;
+
+ // evaluate local function
+ RangeType localFunctionEvaluated(0.0);
+ localFunction.evaluate(localPoint, localFunctionEvaluated);
+
+ // 1.0 * v(x)
+ ret = 1.0 * localFunctionEvaluated;
+
+ // return
+ return ret;
+ }
+
+}; // end class ProductOperation
+
+/**
+ * @class AnalyticalDirichletValues
+ * function representing the dirichlet data g for the poisson problem
+ * a \laplace u = f
+ */
+template <class FunctionSpaceImp>
+class AnalyticalDirichletValues
+{
+public:
+ typedef FunctionSpaceImp FunctionSpaceType;
+
+ typedef typename FunctionSpaceType::DomainType DomainType;
+
+ typedef typename FunctionSpaceType::RangeType RangeType;
+
+ typedef typename FunctionSpaceType::DomainFieldType DomainFieldType;
+
+ typedef typename FunctionSpaceType::RangeFieldType RangeFieldType;
+
+ typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
+
+ void evaluate(const DomainType& arg, RangeType& ret) const
+ {
+ ret = 1.0;
+
+ if (!(arg[0] > 0.0)) {
+ // bottom
+ ret = 1.0;
+ } else if (!(arg[0] < 1.0)) {
+ // top
+ ret = 1.0;
+ }
+ if (!(arg[1] > 0.0)) {
+ // left
+ ret = 1.0;
+ } else if (!(arg[1] < 1.0)) {
+ // right
+ ret = 1.0;
+ }
+ }
+
+ void jacobian(const DomainType& arg, JacobianRangeType& ret) const
+ {
+ ret = 0.0;
+ }
+};
+
+
+// disable warnings about problems in dune headers
+#include <dune/fem-tools/header/disablewarnings.hh>
+
+
+int main(int argc, char** argv)
+{
+ try {
+
+ // MPI manager
+ Dune::MPIManager::initialize(argc, argv);
+
+
+ // grid
+ static const unsigned int dimRange = 1;
+
+ typedef Dune::GridSelector::GridType HGridType;
+
+ typedef Dune::AdaptiveLeafGridPart<HGridType> GridPartType;
+
+ Dune::GridPtr<HGridType> gridPtr("macrogrids/unitcube2.dgf");
+
+ GridPartType gridPart(*gridPtr);
+
+
+ // function spaces and functions
+ typedef Dune::FunctionSpace<double, double, HGridType::dimension, dimRange> FunctionSpaceType;
+
+ typedef Dune::Function<double, double> FunctionType;
+
+ typedef Dune::LagrangeDiscreteFunctionSpace<FunctionSpaceType, GridPartType, polOrder> DiscreteH1Type;
+
+ DiscreteH1Type discreteH1(gridPart);
+
+ typedef Dune::AdaptiveDiscreteFunction<DiscreteH1Type> DiscreteFunctionType;
+
+ typedef Constraints::Dirichlet<DiscreteH1Type> DirichletConstraints;
+
+ DirichletConstraints dirichletConstraints(discreteH1);
+
+ typedef Subspace::Linear<DiscreteH1Type, DirichletConstraints> DiscreteH10Type;
+
+ DiscreteH10Type discreteH10(discreteH1, dirichletConstraints);
+
+ typedef AnalyticalDirichletValues<FunctionSpaceType> AnalyticalDirichletValuesType;
+
+ AnalyticalDirichletValuesType analyticalDirichletValues;
+
+ DiscreteFunctionType discreteDirichletValues("dirichlet_values", discreteH1);
+
+ Dune::FemTools::BetterL2Projection::project(analyticalDirichletValues, discreteDirichletValues);
+
+ // FunctionType should be obsolete
+ // It should be either a DoF-Container or an analytical function, and for
+ // transition phase a discrete function.
+ typedef Subspace::Affine<DiscreteH10Type, DiscreteFunctionType> DiscreteH1gType;
+
+ DiscreteH1gType discreteH1g(discreteH10, discreteDirichletValues);
+
+
+ // operator and functional
+ typedef Operator::FiniteElement<DiscreteH1Type, EllipticOperation> FEMellipticOperator;
+
+ EllipticOperation ellipticOperation;
+
+ FEMellipticOperator femEllipticOperator(discreteH1, ellipticOperation);
+
+ typedef Functional::FiniteElement<DiscreteH1Type, ProductOperation> FEMrhsFunctional;
+
+ ProductOperation productOperation;
+
+ FEMrhsFunctional femRhsFunctional(discreteH1, productOperation);
+
+
+ // matrix, rhs and solution storage
+ typedef Dune::FieldMatrix<double, dimRange, dimRange> FieldMatrixType;
+
+ typedef Container::MatrixFactory<Dune::BCRSMatrix<FieldMatrixType>> MatrixFactoryType;
+
+ typedef MatrixFactoryType::AutoPtrType MatrixContainerPtr;
+
+ typedef MatrixFactoryType::ContainerType MatrixContainer;
+
+ typedef Container::VectorFactory<Dune::BlockVector<Dune::FieldVector<double, 1>>> VectorFactoryType;
+
+ typedef VectorFactoryType::AutoPtrType VectorContainerPtr;
+
+ typedef VectorFactoryType::ContainerType VectorContainer;
+
+ // SparsityPattern & pattern = h1.fullSparsityPattern();
+
+ MatrixContainerPtr A = MatrixFactoryType::create(discreteH10);
+ VectorContainerPtr F = VectorFactoryType::create(discreteH10);
+ VectorContainerPtr G = VectorFactoryType::create(discreteH1);
+
+ // MatrixContainer& A = Container::MatrixFactory<MatrixContainer>::createRef( discreteH10 );
+ // VectorContainer& F = Container::VectorFactory<VectorContainer>::createRef( discreteH10 );
+ // VectorContainer& G = Container::VectorFactory<VectorContainer>::createRef( discreteH1 );
+
+ VectorContainerPtr u0 = VectorFactoryType::create(discreteH10);
+
+
+ // assembler
+ typedef Solver::FEMAssembler<MatrixContainer, VectorContainer> Assembler;
+
+ Assembler::assembleMatrix(femEllipticOperator, *A);
+ Assembler::applyMatrixConstraints(discreteH10, *A);
+
+ Assembler::assembleVector(femRhsFunctional, *F);
+ Assembler::applyVectorConstraints(discreteH10, *F);
+
+ // Assembler::assembleVector( femEllipticOperator( discreteDirichletValues ), *G );
+
+
+ // preconditioner and solver
+ typedef Dune::MatrixAdapter<MatrixContainer, VectorContainer, VectorContainer> MatrixAdapterType;
+
+ MatrixAdapterType op(*A);
+
+ typedef Dune::SeqILU0<MatrixContainer, VectorContainer, VectorContainer, 1> SeqILU0Type;
+
+ SeqILU0Type prec(*A, 1.0);
+
+ typedef Dune::CGSolver<VectorContainer> CG;
+
+ CG cg(op, prec, 1e-4, 100, 2);
+
+ Dune::InverseOperatorResult res;
+
+ // u_0 = A^(-1) ( F - G )
+ // *G = G->operator*=( -1.0 );
+ // *F = F->operator+=( *G );
+ cg.apply(*u0, *F, res);
+
+
+ // postprocessing
+ DiscreteFunctionType solution = Dune::FemTools::discreteFunctionFactory<DiscreteFunctionType>(discreteH1, *u0);
+ // solution += discreteDirichletValues;
+ Dune::FemTools::writeDiscreteFunctionToVTK(solution, "solution");
+
+ return 0;
+ } catch (Dune::Exception& e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ } catch (...) {
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ }
+}
diff --git a/examples/finite_element_local_operation/macrogrids/unitcube1.dgf b/examples/finite_element_local_operation/macrogrids/unitcube1.dgf
new file mode 100644
index 000000000..e3d9007ea
--- /dev/null
+++ b/examples/finite_element_local_operation/macrogrids/unitcube1.dgf
@@ -0,0 +1,18 @@
+DGF
+
+Interval
+0 % first corner
+1 % second corner
+4 % 4 cells only
+#
+
+BoundaryDomain
+default 1 % all boundaries have id 1
+#
+
+GridParameter
+name 1d_Unit_Cube
+overlap 1
+#
+
+# unitcube1.dgf
diff --git a/examples/finite_element_local_operation/macrogrids/unitcube2.dgf b/examples/finite_element_local_operation/macrogrids/unitcube2.dgf
new file mode 100644
index 000000000..7a21d41e1
--- /dev/null
+++ b/examples/finite_element_local_operation/macrogrids/unitcube2.dgf
@@ -0,0 +1,18 @@
+DGF
+
+Interval
+0 0 % first corner
+1 1 % second corner
+40 40 % 4 cells in each direction (x, y)
+#
+
+BoundaryDomain
+default 1 % all boundaries have id 1
+#
+
+GridParameter
+name 2d_Unit_Cube
+overlap 1
+#
+
+# unitcube2.dgf
diff --git a/examples/finite_element_local_operation/macrogrids/unitcube3.dgf b/examples/finite_element_local_operation/macrogrids/unitcube3.dgf
new file mode 100644
index 000000000..1753fa226
--- /dev/null
+++ b/examples/finite_element_local_operation/macrogrids/unitcube3.dgf
@@ -0,0 +1,18 @@
+DGF
+
+Interval
+0 0 0 % first corner
+1 1 1 % second corner
+4 4 4 % 4 cells in each direction (x, y, z)
+#
+
+BoundaryDomain
+default 1 % all boundaries have id 1
+#
+
+GridParameter
+name 3d_Unit_Cube
+overlap 1
+#
+
+# unitcube3.dgf
--
GitLab