diff --git a/src/averageinterface.hh b/src/averageinterface.hh
index e82d08005fca0fd42b3ec635e02534fe693d274f..d2f18ffc86f44d5f65f3e88813ef445b5356b2ac 100644
--- a/src/averageinterface.hh
+++ b/src/averageinterface.hh
@@ -6,10 +6,656 @@
#include <dune/ag-common/dgindexset.hh>
#include <dune/ag-common/crossproduct.hh>
+#include <dune/ag-common/surfmassmatrix.hh>
#include "svd.hh"
#include "lapackpp.h"
#undef max
+template <class GridType>
+class PressureAverager : public Ipopt::TNLP
+{
+ typedef double field_type;
+
+ typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > MatrixType;
+ typedef typename MatrixType::row_type RowType;
+
+
+ enum {dim=GridType::dimension};
+
+public:
+ /** \brief Constructor */
+ PressureAverager(const BoundaryPatch<GridType>* patch,
+ Dune::BlockVector<Dune::FieldVector<double,dim> >* result,
+ const Dune::FieldVector<double,dim>& resultantForce,
+ const Dune::FieldVector<double,dim>& resultantTorque,
+ const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >* massMatrix,
+ const Dune::BlockVector<Dune::FieldVector<double,1> >* nodalWeights,
+ const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >* constraintJacobian)
+ : jacobianCutoff_(1e-8), patch_(patch), x_(result),
+ massMatrix_(massMatrix), nodalWeights_(nodalWeights),
+ constraintJacobian_(constraintJacobian),
+ resultantForce_(resultantForce), resultantTorque_(resultantTorque)
+ {
+ patchArea_ = patch->area();
+ }
+
+ /** default destructor */
+ virtual ~PressureAverager() {};
+
+ /**@name Overloaded from TNLP */
+ //@{
+ /** Method to return some info about the nlp */
+ virtual bool get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+ Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style);
+
+ /** Method to return the bounds for my problem */
+ virtual bool get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+ Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u);
+
+ /** Method to return the starting point for the algorithm */
+ virtual bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+ bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+ Ipopt::Index m, bool init_lambda,
+ Ipopt::Number* lambda);
+
+ /** Method to return the objective value */
+ virtual bool eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value);
+
+ /** Method to return the gradient of the objective */
+ virtual bool eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f);
+
+ /** Method to return the constraint residuals */
+ virtual bool eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g);
+
+ /** Method to return:
+ * 1) The structure of the jacobian (if "values" is NULL)
+ * 2) The values of the jacobian (if "values" is not NULL)
+ */
+ virtual bool eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+ Ipopt::Number* values);
+
+ /** Method to return:
+ * 1) The structure of the hessian of the lagrangian (if "values" is NULL)
+ * 2) The values of the hessian of the lagrangian (if "values" is not NULL)
+ */
+ virtual bool eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+ bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+ Ipopt::Index* jCol, Ipopt::Number* values);
+
+ //@}
+
+ /** @name Solution Methods */
+ //@{
+ /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
+ virtual void finalize_solution(Ipopt::SolverReturn status,
+ Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+ Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+ Ipopt::Number obj_value,
+ const Ipopt::IpoptData* ip_data,
+ Ipopt::IpoptCalculatedQuantities* ip_cq);
+ //@}
+
+ // /////////////////////////////////
+ // Data
+ // /////////////////////////////////
+
+ /** \brief All entries in the constraint Jacobian smaller than the value
+ here are removed. This increases stability.
+ */
+ const double jacobianCutoff_;
+
+ const BoundaryPatch<GridType>* patch_;
+
+ double patchArea_;
+
+ const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >* massMatrix_;
+
+ const Dune::BlockVector<Dune::FieldVector<double,1> >* nodalWeights_;
+
+ const Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >* constraintJacobian_;
+
+ Dune::BlockVector<Dune::FieldVector<double,dim> >* x_;
+
+ Dune::FieldVector<double,dim> resultantForce_;
+ Dune::FieldVector<double,dim> resultantTorque_;
+
+private:
+ /**@name Methods to block default compiler methods.
+ */
+ //@{
+ // PressureAverager();
+ PressureAverager(const PressureAverager&);
+ PressureAverager& operator=(const PressureAverager&);
+ //@}
+};
+
+// returns the size of the problem
+template <class GridType>
+bool PressureAverager<GridType>::
+get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+ Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style)
+{
+ // One variable for each vertex on the coupling boundary, and three for the closest constant field
+ n = patch_->numVertices()*dim + dim;
+
+ // prescribed total forces and moments
+ m = 2*dim;
+
+ // number of nonzeroes in the constraint Jacobian
+ // leave out the very small ones, as they create instabilities
+ nnz_jac_g = 0;
+ for (int i=0; i<m; i++) {
+
+ const RowType& jacobianRow = (*constraintJacobian_)[i];
+
+ for (typename RowType::ConstIterator cIt = jacobianRow.begin(); cIt!=jacobianRow.end(); ++cIt)
+ if ( (*cIt)[0][0] > jacobianCutoff_ )
+ nnz_jac_g++;
+
+ }
+
+ // We only need the lower left corner of the Hessian (since it is symmetric)
+ if (!massMatrix_)
+ DUNE_THROW(SolverError, "No mass matrix has been supplied!");
+
+ nnz_h_lag = 0;
+
+ // use the C style indexing (0-based)
+ index_style = Ipopt::TNLP::C_STYLE;
+
+ return true;
+}
+
+
+// returns the variable bounds
+template <class GridType>
+bool PressureAverager<GridType>::
+get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+ Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u)
+{
+ // here, the n and m we gave IPOPT in get_nlp_info are passed back to us.
+ // If desired, we could assert to make sure they are what we think they are.
+ //assert(n == x_->dim());
+ //assert(m == 0);
+
+ // Be on the safe side: unset all variable bounds
+ for (size_t i=0; i<n; i++) {
+ x_l[i] = -std::numeric_limits<double>::max();
+ x_u[i] = std::numeric_limits<double>::max();
+ }
+
+ for (int i=0; i<dim; i++) {
+ g_l[i] = g_u[i] = resultantForce_[i];
+ g_l[i+dim] = g_u[i+dim] = resultantTorque_[i];
+ }
+
+ return true;
+}
+
+// returns the initial point for the problem
+template <class GridType>
+bool PressureAverager<GridType>::
+get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+ bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+ Ipopt::Index m, bool init_lambda, Ipopt::Number* lambda)
+{
+ // Here, we assume we only have starting values for x, if you code
+ // your own NLP, you can provide starting values for the dual variables
+ // if you wish
+ assert(init_x == true);
+ assert(init_z == false);
+ assert(init_lambda == false);
+
+ // initialize to the given starting point
+ for (int i=0; i<n; i++)
+ x[i] = 0;
+
+ return true;
+}
+
+// returns the value of the objective function
+template <class GridType>
+bool PressureAverager<GridType>::
+eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value)
+{
+// std::cout << "x:" << std::endl;
+// for (int i=0; i<n; i++)
+// std::cout << x[i] << std::endl;
+
+ // Init return value
+ obj_value = 0;
+
+ ////////////////////////////////////
+ // Compute x^T*A*x
+ ////////////////////////////////////
+
+ for (int rowIdx=0; rowIdx<massMatrix_->N(); rowIdx++) {
+
+ const typename MatrixType::row_type& row = (*massMatrix_)[rowIdx];
+
+ typename MatrixType::row_type::ConstIterator cIt = row.begin();
+ typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt)
+ for (int i=0; i<dim; i++)
+ obj_value += x[dim*rowIdx+i] * x[dim*cIt.index()+i] * (*cIt)[0][0];
+
+ }
+
+ // -b(x)
+ for (int i=0; i<nodalWeights_->size(); i++)
+ for (int j=0; j<dim; j++)
+ obj_value -= 2 * x[n-dim + j] * x[i*dim+j] * (*nodalWeights_)[i];
+
+ // += c^2 * \int 1 ds
+ for (int i=0; i<dim; i++)
+ obj_value += patchArea_ * (x[n-dim + i] * x[n-dim + i]);
+
+ //std::cout << "IPOPT Energy: " << obj_value << std::endl;
+ //exit(0);
+
+ return true;
+}
+
+// return the gradient of the objective function grad_{x} f(x)
+template <class GridType>
+bool PressureAverager<GridType>::
+eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f)
+{
+ //std::cout << "### eval_grad_f ###" << std::endl;
+
+ // \nabla J = A(x,.) - b(x)
+ for (int i=0; i<n; i++)
+ grad_f[i] = 0;
+
+ printmatrix(std::cout, *massMatrix_, "mass", "--");
+
+ for (int i=0; i<massMatrix_->N(); i++) {
+
+ const typename MatrixType::row_type& row = (*massMatrix_)[i];
+
+ typename MatrixType::row_type::ConstIterator cIt = row.begin();
+ typename MatrixType::row_type::ConstIterator cEndIt = row.end();
+
+ for (; cIt!=cEndIt; ++cIt)
+ for (int j=0; j<dim; j++)
+ grad_f[i*dim+j] += 2 * (*cIt)[0][0] * x[cIt.index()*dim+j];
+
+ for (int j=0; j<dim; j++)
+ grad_f[i*dim+j] -= 2 * x[n-dim+j] * (*nodalWeights_)[i];
+
+ }
+
+ for (int i=0; i<dim; i++) {
+
+ for (int j=0; j<nodalWeights_->size(); j++)
+ grad_f[n-dim+i] -= 2* (*nodalWeights_)[j]*x[j*dim+i];
+
+ grad_f[n-dim+i] += 2*x[n-dim+i]*patchArea_;
+
+ }
+
+ for (int i=0; i<n; i++) {
+ std::cout << "x = " << x[i] << std::endl;
+ std::cout << "grad = " << grad_f[i] << std::endl;
+ }
+
+ return true;
+}
+
+// return the value of the constraints: g(x)
+template <class GridType>
+bool PressureAverager<GridType>::
+eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+ for (int i=0; i<m; i++) {
+
+ // init
+ g[i] = 0;
+
+ const RowType& jacobianRow = (*constraintJacobian_)[i];
+
+ for (typename RowType::ConstIterator cIt = jacobianRow.begin(); cIt!=jacobianRow.end(); ++cIt)
+ if ( (*cIt)[0][0] > jacobianCutoff_ )
+ g[i] += (*cIt)[0][0] * x[cIt.index()];
+
+ }
+
+ return true;
+}
+
+// return the structure or values of the jacobian
+template <class GridType>
+bool PressureAverager<GridType>::
+eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+ Ipopt::Number* values)
+{
+ int idx = 0;
+
+ if (values==NULL) {
+
+ for (int i=0; i<m; i++) {
+
+ const RowType& jacobianRow = (*constraintJacobian_)[i];
+
+ for (typename RowType::ConstIterator cIt = jacobianRow.begin(); cIt!=jacobianRow.end(); ++cIt) {
+ if ( (*cIt)[0][0] > jacobianCutoff_ ) {
+ iRow[idx] = i;
+ jCol[idx] = cIt.index();
+ idx++;
+ }
+ }
+
+ }
+
+ } else {
+
+ for (int i=0; i<m; i++) {
+
+ const RowType& jacobianRow = (*constraintJacobian_)[i];
+
+ for (typename RowType::ConstIterator cIt = jacobianRow.begin(); cIt!=jacobianRow.end(); ++cIt)
+ if ( (*cIt)[0][0] > jacobianCutoff_ )
+ values[idx++] = (*cIt)[0][0];
+
+ }
+
+
+
+ }
+
+ return true;
+}
+
+//return the structure or values of the hessian
+template <class GridType>
+bool PressureAverager<GridType>::
+eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+ Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+ bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+ Ipopt::Index* jCol, Ipopt::Number* values)
+{
+ // We are using a quasi-Hessian approximation
+ return false;
+}
+
+template <class GridType>
+void PressureAverager<GridType>::
+finalize_solution(Ipopt::SolverReturn status,
+ Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+ Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+ Ipopt::Number obj_value,
+ const Ipopt::IpoptData* ip_data,
+ Ipopt::IpoptCalculatedQuantities* ip_cq)
+{
+ x_->resize(patch_->numVertices());
+
+ for (int i=0; i<x_->size(); i++)
+ for (int j=0; j<dim; j++)
+ (*x_)[i][j] = x[i*dim+j];
+
+ std::cout << "Closest constant: ";// << x[n-dim] << " " << x[n-dim+1] << " " << [x-dim+2] << std::endl;
+
+}
+
+
+
+// Given a resultant force and torque (from a rod problem), this method computes the corresponding
+// Neumann data for a 3d elasticity problem.
+template <class GridType>
+void computeAveragePressureIPOpt(const Dune::FieldVector<double,GridType::dimension>& resultantForce,
+ const Dune::FieldVector<double,GridType::dimension>& resultantTorque,
+ const BoundaryPatch<GridType>& interface,
+ const Configuration& crossSection,
+ Dune::BlockVector<Dune::FieldVector<double, GridType::dimension> >& pressure)
+{
+ const GridType& grid = interface.getGrid();
+ const int level = interface.level();
+ const typename GridType::Traits::LevelIndexSet& indexSet = grid.levelIndexSet(level);
+ const int dim = GridType::dimension;
+ typedef typename GridType::ctype ctype;
+ typedef double field_type;
+
+ typedef typename GridType::template Codim<dim>::LevelIterator VertexIterator;
+
+ // Create the matrix of constraints
+ Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> > matrix(2*dim, dim*interface.numVertices(),
+ Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >::random);
+
+ for (int i=0; i<dim; i++) {
+ matrix.setrowsize(i, interface.numVertices());
+ matrix.setrowsize(i+dim, dim*interface.numVertices());
+ }
+
+ matrix.endrowsizes();
+
+ for (int i=0; i<dim; i++)
+ for (int j=0; j<interface.numVertices(); j++)
+ matrix.addindex(i, dim*j+i);
+
+ for (int i=0; i<dim; i++)
+ for (int j=0; j<dim*interface.numVertices(); j++)
+ matrix.addindex(i+dim, j);
+
+ matrix.endindices();
+
+ matrix = 0;
+
+ // Create the surface mass matrix
+ Dune::BCRSMatrix<Dune::FieldMatrix<field_type,1,1> > massMatrix;
+ assembleSurfaceMassMatrix<GridType,1>(interface, massMatrix);
+
+ // Make global-to-local array
+ std::vector<int> globalToLocal;
+ interface.makeGlobalToLocal(globalToLocal);
+
+ // Make array of nodal weights
+ Dune::BlockVector<Dune::FieldVector<double,1> > nodalWeights(interface.numVertices());
+ nodalWeights = 0;
+
+ typename GridType::template Codim<0>::LevelIterator eIt = indexSet.template begin<0,Dune::All_Partition>();
+ typename GridType::template Codim<0>::LevelIterator eEndIt = indexSet.template end<0,Dune::All_Partition>();
+
+ for (; eIt!=eEndIt; ++eIt) {
+
+ typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nIt = eIt->ilevelbegin();
+ typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nEndIt = eIt->ilevelend();
+
+ for (; nIt!=nEndIt; ++nIt) {
+
+ if (!interface.contains(*eIt,nIt))
+ continue;
+
+ const Dune::LagrangeShapeFunctionSet<ctype, field_type, dim-1>& baseSet
+ = Dune::LagrangeShapeFunctions<ctype, field_type, dim-1>::general(nIt.intersectionGlobal().type(),1);
+
+ const Dune::ReferenceElement<double,dim>& refElement = Dune::ReferenceElements<double, dim>::general(eIt->type());
+
+ // four rows because a face may have no more than four vertices
+ Dune::FieldVector<double,4> mu(0);
+ Dune::FieldVector<double,3> mu_tilde[4][3];
+
+ for (int i=0; i<4; i++)
+ for (int j=0; j<3; j++)
+ mu_tilde[i][j] = 0;
+
+ for (int i=0; i<nIt.intersectionGlobal().corners(); i++) {
+
+ const Dune::QuadratureRule<double, dim-1>& quad
+ = Dune::QuadratureRules<double, dim-1>::rule(nIt.intersectionGlobal().type(), dim-1);
+
+ for (size_t qp=0; qp<quad.size(); qp++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<double,dim-1>& quadPos = quad[qp].position();
+
+ const double integrationElement = nIt.intersectionGlobal().integrationElement(quadPos);
+
+ // \mu_i = \int_t \varphi_i \ds
+ mu[i] += quad[qp].weight() * integrationElement * baseSet[i].evaluateFunction(0,quadPos);
+
+ // \tilde{\mu}_i^j = \int_t \varphi_i \times (x - x_0) \ds
+ Dune::FieldVector<double,dim> worldPos = nIt.intersectionGlobal().global(quadPos);
+
+ for (int j=0; j<dim; j++) {
+
+ // Vector-valued basis function
+ Dune::FieldVector<double,dim> phi_i(0);
+ phi_i[j] = baseSet[i].evaluateFunction(0,quadPos);
+
+ mu_tilde[i][j].axpy(quad[qp].weight() * integrationElement,
+ crossProduct(worldPos-crossSection.r, phi_i));
+
+ }
+
+ }
+
+ }
+
+ // Set up matrix
+ for (int i=0; i<baseSet.size(); i++) {
+
+ int faceIdxi = refElement.subEntity(nIt.numberInSelf(), 1, i, dim);
+ int subIndex = globalToLocal[indexSet.template subIndex<dim>(*eIt, faceIdxi)];
+
+ nodalWeights[subIndex] += mu[i];
+ for (int j=0; j<dim; j++)
+ matrix[j][subIndex*dim+j] += mu[i];
+
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ matrix[dim+k][dim*subIndex+j] += mu_tilde[i][j][k];
+
+ }
+
+ }
+
+ }
+
+ //printmatrix(std::cout, matrix, "jacobian", "--");
+ //printmatrix(std::cout, massMatrix, "mass", "--");
+
+ // /////////////////////////////////////////////////////////////////////////////////////
+ // Set up and start the interior-point solver
+ // /////////////////////////////////////////////////////////////////////////////////////
+
+ // Create a new instance of IpoptApplication
+ Ipopt::SmartPtr<Ipopt::IpoptApplication> app = new Ipopt::IpoptApplication();
+
+ // Change some options
+ app->Options()->SetNumericValue("tol", 1e-8);
+ app->Options()->SetIntegerValue("max_iter", 20);
+ app->Options()->SetStringValue("mu_strategy", "adaptive");
+ app->Options()->SetStringValue("output_file", "ipopt.out");
+ app->Options()->SetStringValue("hessian_approximation", "limited-memory");
+
+ // Intialize the IpoptApplication and process the options
+ Ipopt::ApplicationReturnStatus status;
+ status = app->Initialize();
+ if (status != Ipopt::Solve_Succeeded)
+ DUNE_THROW(SolverError, "Error during IPOpt initialization!");
+
+ // Ask Ipopt to solve the problem
+ Dune::BlockVector<Dune::FieldVector<double,dim> > localPressure;
+ Ipopt::SmartPtr<Ipopt::TNLP> defectSolverSmart = new PressureAverager<GridType>(&interface,
+ &localPressure,
+ resultantForce,
+ resultantTorque,
+ &massMatrix,
+ &nodalWeights,
+ &matrix);
+ status = app->OptimizeTNLP(defectSolverSmart);
+
+ if (status != Ipopt::Solve_Succeeded)
+ DUNE_THROW(SolverError, "Solving the defect problem failed!");
+
+ // //////////////////////////////////////////////////////////////////////////////
+ // Get result
+ // //////////////////////////////////////////////////////////////////////////////
+
+ // set up output array
+ pressure.resize(indexSet.size(dim));
+ pressure = 0;
+
+ for (size_t i=0; i<globalToLocal.size(); i++)
+ if (globalToLocal[i]>=0)
+ pressure[i] = localPressure[globalToLocal[i]];
+
+ // /////////////////////////////////////////////////////////////////////////////////////
+ // Compute the overall force and torque to see whether the preceding code is correct
+ // /////////////////////////////////////////////////////////////////////////////////////
+#if 1
+ Dune::FieldVector<double,3> outputForce(0), outputTorque(0);
+
+ eIt = indexSet.template begin<0,Dune::All_Partition>();
+ eEndIt = indexSet.template end<0,Dune::All_Partition>();
+
+ for (; eIt!=eEndIt; ++eIt) {
+
+ typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nIt = eIt->ilevelbegin();
+ typename GridType::template Codim<0>::Entity::LevelIntersectionIterator nEndIt = eIt->ilevelend();
+
+ for (; nIt!=nEndIt; ++nIt) {
+
+ if (!interface.contains(*eIt,nIt))
+ continue;
+
+ const Dune::LagrangeShapeFunctionSet<double, double, dim-1>& baseSet
+ = Dune::LagrangeShapeFunctions<double, double, dim-1>::general(nIt.intersectionGlobal().type(),1);
+
+ const Dune::QuadratureRule<double, dim-1>& quad
+ = Dune::QuadratureRules<double, dim-1>::rule(nIt.intersectionGlobal().type(), dim-1);
+
+ const Dune::ReferenceElement<double,dim>& refElement = Dune::ReferenceElements<double, dim>::general(eIt->type());
+
+ for (size_t qp=0; qp<quad.size(); qp++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<double,dim-1>& quadPos = quad[qp].position();
+
+ const double integrationElement = nIt.intersectionGlobal().integrationElement(quadPos);
+
+ // Evaluate function
+ Dune::FieldVector<double,dim> localPressure(0);
+
+ for (size_t i=0; i<baseSet.size(); i++) {
+
+ int faceIdxi = refElement.subEntity(nIt.numberInSelf(), 1, i, dim);
+ int subIndex = indexSet.template subIndex<dim>(*eIt, faceIdxi);
+
+ localPressure.axpy(baseSet[i].evaluateFunction(0,quadPos),
+ pressure[subIndex]);
+
+ }
+
+ // Sum up the total force
+ outputForce.axpy(quad[qp].weight()*integrationElement, localPressure);
+
+ // Sum up the total torque \int (x - x_0) \times f dx
+ Dune::FieldVector<double,dim> worldPos = nIt.intersectionGlobal().global(quadPos);
+ outputTorque.axpy(quad[qp].weight()*integrationElement,
+ crossProduct(worldPos - crossSection.r, localPressure));
+
+ }
+
+ }
+
+ }
+
+ outputForce -= resultantForce;
+ outputTorque -= resultantTorque;
+ assert( outputForce.infinity_norm() < 1e-6 );
+ assert( outputTorque.infinity_norm() < 1e-6 );
+// std::cout << "Output force: " << outputForce << std::endl;
+// std::cout << "Output torque: " << outputTorque << " " << resultantTorque[0]/outputTorque[0] << std::endl;
+#endif
+
+}
+
+
// Given a resultant force and torque (from a rod problem), this method computes the corresponding
// Neumann data for a 3d elasticity problem.
template <class GridType>
@@ -145,6 +791,12 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
b(i+3) = resultantTorque[i] * segmentArea;
}
+ for (int i=0; i<6; i++) {
+ for (int j=0; j<3*baseSet.size(); j++)
+ std::cout << matrix(i,j) << " ";
+ std::cout << std::endl;
+ }
+
LaLinearSolve(matrix, u, b);
#endif
// std::cout << b << std::endl;
@@ -224,6 +876,8 @@ void computeAveragePressure(const Dune::FieldVector<double,GridType::dimension>&
}
+
+
template <class GridType>
void averageSurfaceDGFunction(const GridType& grid,
const Dune::BlockVector<Dune::FieldVector<double,GridType::dimension> >& dgFunction,