// compute inner part of (global) stiffnes matrix
  // 
  template <class PDE, class INT, class MATRIX, class MESH>
  void element_compute_stima(MATRIX& A, 
                             const MESH& M, 
                             const PDE& P,
                             const INT& integration_order,
                             const typename MESH::size_type level,
                             std::pair<typename MESH::size_type, typename MESH::size_type> element_range )
  {
    // loop over all elements and sum up local stiffness matrices
    const typename MESH::size_type              MM = MESH::nodes_per_element;

    typename MESH::local_stiffness_matrix_type  stima(MM,MM); 

#ifdef USE_INSERT
    typename MESH::size_type sortflag = (M.element_count()/4)+1;
#endif

#ifndef USEPARALLELMATRIX
#pragma omp parallel for default(shared) firstprivate(stima) schedule(static)
#endif
    for (typename MESH::size_type ei=element_range.first; ei < element_range.second; ++ei) {
      if (M.element_hierarchy(ei).level > level) continue;
      if (M.element_hierarchy(ei).child[0] != MESH::NOCHILD 
          && M.element_hierarchy(M.element_hierarchy(ei).child[0]).level <= level) continue;
      stima.clear();
      compute_local_stima( stima, 
                           M, 
                           ei, 
                           typename MESH::element_type(), 
                           P,
                           integration_order );
      // update global stima
      for (typename MESH::size_type j=0; j<MM; ++j) {
        const typename MESH::size_type nj = M.element(ei)(j);
#ifndef USEPARALLELMATRIX
#pragma omp critical (element_compute_stima_update_global)
#endif
        for (typename MESH::size_type i=0; i<MM; ++i) {
#ifndef USE_INSERT
          A(M.element(ei)(i), nj) += stima(i,j);
#else
          A.append_element(M.element(ei)(i), nj, stima(i,j));
#endif
        }
      }
#ifdef USE_INSERT
      if ( --sortflag == 0) { 
#ifndef USEPARALLELMATRIX
#pragma omp critical (element_compute_stima_update_global)
#endif
        {
          A.sort(); 
          sortflag = (M.element_count()/4)+1; 
        }
      }
#endif
    }
    
    return;
  }