diff --git a/ugbase/lib_disc/operator/linear_operator/subspace_correction/sequential_subspace_correction.h b/ugbase/lib_disc/operator/linear_operator/subspace_correction/sequential_subspace_correction.h index 5b7387fe6..2e2bb9e0b 100644 --- a/ugbase/lib_disc/operator/linear_operator/subspace_correction/sequential_subspace_correction.h +++ b/ugbase/lib_disc/operator/linear_operator/subspace_correction/sequential_subspace_correction.h @@ -54,7 +54,7 @@ namespace ug{ //! Abstract definition for subspace V_k -template +template class ILocalSubspace { public: @@ -67,28 +67,32 @@ class ILocalSubspace /// Matrix type typedef typename TAlgebra::matrix_type matrix_type; + /// Grid function type + typedef GridFunction grid_function_type; typedef DenseVector< VariableArray1 > vector_type_local; typedef DenseMatrix< VariableArray2 > matrix_type_local; /// virtual DTOR - virtual ~ILocalSubspace(){} + virtual ~ILocalSubspace(){} /// Called once. - virtual bool preprocess(const vector_type &c) {return true;} + virtual bool preprocess(const vector_type &c) + { return true; } - /// Extract local data (based on group obj) - virtual void init(TGroupObj*, const vector_type &) = 0; - /// Extract matrix (on local index set) + /// Extract local data (based on group obj). + virtual void init(TObject*, const vector_type &) = 0; + + /// Extract matrix (on local index set). virtual void extract_matrix(const matrix_type &A) = 0; - /// Extract rhs (on local index set) for parallel subspace correction + /// Extract rhs (on local index set) for parallel subspace correction. virtual void extract_rhs(const vector_type &d) = 0; - /// Extract rhs (on local index set) for sequential subspace correction + /// Extract rhs (on local index set) for sequential subspace correction. virtual void extract_rhs(const vector_type &d, const matrix_type &A, const vector_type &c) = 0; - /// u = u + omega*ck + /// u = u + omega*ck. virtual void update_solution(vector_type &u, double omega=1.0) = 0; virtual size_t size() { return 0; } @@ -96,9 +100,9 @@ class ILocalSubspace -//! Abstract definition for subspace V_k -template -class LocalIndexSubspace : public ILocalSubspace +//! Concrete definition of subspace V_k (based on size_t) +template +class LocalIndexSubspace : public ILocalSubspace { public: /// Algebra type @@ -116,13 +120,18 @@ class LocalIndexSubspace : public ILocalSubspace typedef DenseMatrix< VariableArray2 > matrix_type_local; /// virtual DTOR - virtual ~LocalIndexSubspace(){} + virtual ~LocalIndexSubspace(){} /// Called once. - virtual bool preprocess(const vector_type &c) {return true;} - + /*virtual bool preprocess(const vector_type &c, obj_iterator_type &objBegin, obj_iterator_type &objEnd) + { + objBegin = c.template begin(); + objEnd = c.template end(); + return true; + } +*/ /// Extract local data (based on group obj) - virtual void init(TGroupObj*, const vector_type &) = 0; + virtual void init(TObject*, const vector_type &) = 0; /// Extract matrix (on local index set) virtual void extract_matrix(const matrix_type &A) @@ -218,9 +227,9 @@ class LocalIndexSubspace : public ILocalSubspace }; -//! Abstract definition for subspace V_k -template -class LocalDoFSubspace : public ILocalSubspace +//! Abstract definition for subspace V_k (based on DoFIndex) +template +class LocalDoFSubspace : public ILocalSubspace { public: /// Algebra type @@ -238,13 +247,13 @@ class LocalDoFSubspace : public ILocalSubspace typedef DenseMatrix< VariableArray2 > matrix_type_local; /// virtual DTOR - virtual ~LocalDoFSubspace(){} + virtual ~LocalDoFSubspace(){} - /// Called once. - virtual bool preprocess(const vector_type &c) {return true;} + bool check(void *obj) const + { return (dynamic_cast (obj) != NULL); } /// Extract local data (based on group obj) - virtual void init(TGroupObj*, const vector_type &) = 0; + virtual void init(TObject*, const vector_type &) = 0; /// Extract matrix (on local index set) virtual void extract_matrix(const matrix_type &A) @@ -336,8 +345,11 @@ class LocalDoFSubspace : public ILocalSubspace //! Collects indices on all elements with v \in Vtx(elem) template -class VertexBasedSubspace : public LocalIndexSubspace +class VertexBasedSubspace : public LocalIndexSubspace { +protected: + typedef Vertex TObject; + public: /// Algebra type typedef TAlgebra algebra_type; @@ -357,9 +369,15 @@ class VertexBasedSubspace : public LocalIndexSubspace /// virtual DTOR virtual ~VertexBasedSubspace(){} + bool check(void *obj) const + { return true; /*return (dynamic_cast (obj) != NULL);*/ } + /// Extract indices for local DoFs. - void init(Vertex *groupObj, const vector_type &cvec) + void init(void *obj, const vector_type &cvec) { + UG_ASSERT(check(obj), "HUHH: Expecting vertex!"); + Vertex *groupObj = reinterpret_cast (obj); + // We will modify index list of base class. std::vector &vInd = base_type::m_vInd; vInd.clear(); @@ -402,6 +420,8 @@ class VertexCenteredVankaSubspace : public LocalDoFSubspace grid_function_type; + /// Base type typedef LocalDoFSubspace base_type; @@ -413,17 +433,14 @@ class VertexCenteredVankaSubspace : public LocalDoFSubspace(){} /// Extracts function IDs. - bool preprocess(const vector_type &cvec) +protected: + void extract_ids(const grid_function_type *c) { - typedef GridFunction TGridFunction; - const TGridFunction *c = dynamic_cast (&cvec); // Need a grid function here! - UG_COND_THROW(c==NULL, "Requiring a grid function here!"); - ConstSmartPtr ddinfo = c->approx_space()->dof_distribution_info(); UG_COND_THROW(ddinfo.invalid(), "Requiring valid ddinfo!"); - // Vertex functions + // Vertex functions. m_vVtxFct.reserve(m_vVtxCmp.size()); m_vVtxFct.clear(); for(size_t i = 0; i < m_vVtxCmp.size(); ++i) @@ -434,11 +451,21 @@ class VertexCenteredVankaSubspace : public LocalDoFSubspacefct_id_by_name(m_vElemCmp[i].c_str())); - + } +public: + /// OVERRIDE + bool preprocess(const vector_type &cvec) + { + typedef GridFunction TGridFunction; + const TGridFunction *c = dynamic_cast (&cvec); // Need a grid function here! + + UG_COND_THROW(c==NULL, "Requiring a grid function here!"); + extract_ids(c); + return true; } - /// Extract indices for local DoFs. + /// Extract indices for single vertex. void init(Vertex *groupObj, const vector_type &cvec) { // We will modify index list of base class. @@ -545,7 +572,7 @@ class LocalSubspaceCorrection : */ - +//! Abstract loop template void ParallelSubspaceCorrectionLoop(const typename TAlgebra::matrix_type& A, GridFunction& c, @@ -569,22 +596,73 @@ void ParallelSubspaceCorrectionLoop(const typename TAlgebra::matrix_type& A, } } -template +//! Abstract loop over TGroupObj (forward) +template void SequentialSubspaceCorrectionLoop(const typename TAlgebra::matrix_type& A, GridFunction& c, const typename TAlgebra::vector_type& d, number omega_relax, - ILocalSubspace &subspace, + ILocalSubspace &subspace, + TGroupObjIter objIterBegin, TGroupObjIter objIterEnd) +{ + // Loop over all grouping objects. + size_t count=0; + for(TGroupObjIter iter = objIterBegin; iter != objIterEnd; ++iter) + { + // Apply subspace correction (w.r.t. obj.) + TGroupObj* groupObj = *iter; + if (groupObj == NULL) + { + UG_LOG("WARNING: Found empty object!: "<< count); + } + else + { + subspace.init(groupObj, c); + subspace.extract_matrix(A); + subspace.extract_rhs(d,A,c); // w/ updates for c => sequential + subspace.update_solution(c, omega_relax); + } + count++; + } +} + +//! Abstract loop (backward) +template +void SequentialSubspaceCorrectionLoopBackward(const typename TAlgebra::matrix_type& A, + GridFunction& c, + const typename TAlgebra::vector_type& d, + number omega_relax, + ILocalSubspace &subspace, typename GridFunction::template traits::const_iterator objIterBegin, typename GridFunction::template traits::const_iterator objIterEnd ) { // Loop over all grouping objects. typedef typename GridFunction::template traits::const_iterator GroupObjIter; + + int ncount = 0; + for(GroupObjIter iter = objIterBegin; iter != objIterEnd; ++iter) + { ncount++; } + + // std::cerr << "Found " << ncount << "objects!" << std::endl; + + std::vector objects(ncount); for(GroupObjIter iter = objIterBegin; iter != objIterEnd; ++iter) + { + // std::cerr << "Found1 " << *iter << std::endl; + objects.push_back(*iter); + } + + // std::cerr << "Created vector !" << std::endl; + + for(typename std::vector::reverse_iterator riter = objects.rbegin(); riter != objects.rend(); ++riter) { // Apply subspace correction (w.r.t. obj.) - TGroupObj* groupObj = *iter; + TGroupObj* groupObj = *riter; + if (groupObj == NULL) continue; + + //std::cerr << "Found2 " << groupObj << std::endl; + subspace.init(groupObj, c); subspace.extract_matrix(A); subspace.extract_rhs(d,A,c); // w/ updates for c => sequential @@ -592,6 +670,85 @@ void SequentialSubspaceCorrectionLoop(const typename TAlgebra::matrix_type& A, } } +template +class customLexLess{ + public: + customLexLess(const typename TDomain::position_accessor_type& aaPos): _aaPos(aaPos){} + bool operator()(Vertex* a, Vertex *b) const + { return _aaPos[a] < _aaPos[b]; } + protected: + const typename TDomain::position_accessor_type& _aaPos; +}; + + +/// Abstract +template +class ISpaceDecomposition { + +public: + /// Algebra type + typedef TAlgebra algebra_type; + + /// Vector type + typedef typename TAlgebra::vector_type vector_type; + + /// Matrix type + typedef typename TAlgebra::matrix_type matrix_type; + + /// Grid function type + typedef GridFunction grid_function_type; + + + virtual ~ISpaceDecomposition(){} + + /// Initialize decomposition + virtual void init(const vector_type &x) = 0; + + /// World Dimension + static const int dim = TDomain::dim; + +}; + + +template +class FullVertexCover: public ISpaceDecomposition +{ +public: + + typedef ISpaceDecomposition base_type; + + typedef GridFunction TGridFunction; + + static const int dim = TDomain::dim; + + typedef typename TGridFunction::element_type TElement; + + FullVertexCover() : m_pSol(NULL) {} + + /// virtual DTOR + virtual ~FullVertexCover(){} + + virtual void init(const typename base_type::vector_type &x) + { + TGridFunction* m_pSol = dynamic_cast(&x); + UG_ASSERT(m_pSol !=NULL, "Error: Cast failed!"); + } + + template + typename TGridFunction::template traits::const_iterator + begin() const + { return m_pSol->template begin();} + + template + typename TGridFunction::template traits::const_iterator + end() const + { return m_pSol->template end(); } +protected: + TGridFunction* m_pSol; + +}; + + /// Sequential subspace correction preconditioner template class SequentialSubspaceCorrection : public IPreconditioner @@ -609,6 +766,9 @@ class SequentialSubspaceCorrection : public IPreconditioner /// Matrix Operator type typedef typename IPreconditioner::matrix_operator_type matrix_operator_type; + /// Subspace for vertices + typedef ILocalSubspace TVertexSubspace; + /// Base type typedef IPreconditioner base_type; @@ -654,7 +814,7 @@ class SequentialSubspaceCorrection : public IPreconditioner void set_type(const std::string& type){ m_type=type; } /// set subspace - void set_vertex_subspace(SmartPtr > spVertexSubspace) + void set_vertex_subspace(SmartPtr spVertexSubspace) { m_spVertexSubspace = spVertexSubspace; } @@ -700,12 +860,14 @@ class SequentialSubspaceCorrection : public IPreconditioner { PROFILE_BEGIN_GROUP(SSC_step, "algebra ssc"); - GridFunction* pC = dynamic_cast*>(&c); + typedef GridFunction TGridFunction; + TGridFunction* pC = dynamic_cast(&c); + UG_COND_THROW(pC == NULL, "SequentialSubspaceCorrection expects correction to be a GridFunction."); //typedef typename GridFunction::element_type Element; //typedef typename GridFunction::side_type Side; - m_spVertexSubspace->preprocess(*pC); + // Set all vector entries to zero. pC->set(0.0); @@ -720,10 +882,53 @@ class SequentialSubspaceCorrection : public IPreconditioner #endif { matrix_type &A=*pOp; - if (m_type == "vertex") SequentialSubspaceCorrectionLoop(A, *pC, d, m_relax, - *m_spVertexSubspace, pC->template begin(), pC->template end()); + if (m_type == "vertex") + { + // Prepare + typedef typename GridFunction::template traits::const_iterator TVertexIter; + TVertexIter objIterBegin = pC->template begin(); + TVertexIter objIterEnd = pC->template end(); + + m_spVertexSubspace->preprocess(*pC); + + if (true) + { + int ncount = 0; + for(TVertexIter iter = objIterBegin; iter != objIterEnd; ++iter) + { ncount++; } + + // std::cerr << "Found " << ncount << "objects!" << std::endl; + + // Create bidirectional container. + typedef std::vector TObjVector; + TObjVector objects(ncount); + objects.clear(); + for(TVertexIter iter = objIterBegin; iter != objIterEnd; ++iter) + { + // std::cerr << "Found1 " << *iter << std::endl; + objects.push_back(*iter); + } + + // Forward loop. + SequentialSubspaceCorrectionLoop + (A, *pC, d, m_relax, *m_spVertexSubspace, objects.begin(), objects.end()); + + // Backward loop. + SequentialSubspaceCorrectionLoop + (A, *pC, d, m_relax, *m_spVertexSubspace, objects.rbegin(), objects.rend()); + + } // true + + + } else if(m_type == "element") { + + /* typedef typename GridFunction::template traits::const_iterator TElemIter; + SequentialSubspaceCorrectionLoop + (A, *pC, d, m_relax, *m_spElemSubspace, pC->template begin(), pC->template end()); +*/ + } else UG_THROW("SequentialSubspaceCorrectionStep: wrong patch type '"<set_storage_type(PST_CONSISTENT); #endif @@ -737,14 +942,15 @@ class SequentialSubspaceCorrection : public IPreconditioner virtual bool postprocess() {return true;} - - +public: + typedef typename GridFunction::element_type Element; protected: number m_relax; std::string m_type; - SmartPtr > m_spVertexSubspace; + SmartPtr m_spVertexSubspace; + #ifdef UG_PARALLEL /// matrix with overlap