[dune-fem] Possible bug when defining a linear operator between spaces whihc live on different grid

Agnese, Marco m.agnese13 at imperial.ac.uk
Thu Feb 26 19:34:00 CET 2015


Hi Dune,
I don't know if I am doing something wrong or if there is a bug in Dune::Fem::Operator. 

I have created an operator derived from Dune::Fem::Operator<DomainFunctionType,RangeFunctionType>.

DomainFunctionType is a P1 continuos lagrangian function from R^2 to R^1 defined over interfacegrid (griddim=1).
RangeFunctionType is a P2 continuos lagrangian function from R^2 to R^2 defined over bulkgrid (griddim=2).

I have implemented the following method inside my operator.

  void debug()
  {
    Stencil<DomainSpaceType,RangeSpaceType> stencil(domainspace_,rangespace_);
    op_.reserve(stencil);
    op_.clear();

    // perform an interface walkthrough and assemble the global matrix
    typedef typename DomainSpaceType::IteratorType InterfaceIteratorType;
    for(InterfaceIteratorType interfaceIt=domainspace_.begin();interfaceIt!=domainspace_.end();++interfaceIt)
    {
      // extract the corresponding bulk entity iterator
      const BulkIteratorType& bulkIt(mapper_.interface2bulk(*interfaceIt));

      // extract local matrix
      typename LinearOperatorType::LocalMatrixType localMatrix(op_.localMatrix(*interfaceIt,*bulkIt));
      const std::size_t columnLocalSize(localMatrix.columns());
      const std::size_t rowLocalSize(localMatrix.rows());
      for(std::size_t i=0;i!=rowLocalSize;++i)
        for(std::size_t j=0;j!=columnLocalSize;++j)
          localMatrix.add(i,j,1.0);
    }
  }

The global matrix has the correct size but all the non-zeros entries are located in the last columns (which is obviously wrong). 
Indeed I should have at least one non-zeros entries in each column. I think that the mapping from the local DOF to global DOF is broken for what concern the column index.
For what concern the rows, it should be correct.

I am using a Dune::Fem::SparseRowLinearOperator while the discrete functions are Dune::Fem::ISTLBlockVectorDiscreteFunction.

I don't know how to fix the problem.

Thank you very much for your help,
best regards,
Marco.



More information about the dune-fem mailing list