| 16 | |
| 17 | template<class SCAL, int DIMA, int DIMB> |
| 18 | shared_ptr<BaseMatrix> ConvertOperator (shared_ptr<FESpace> space_a, shared_ptr<FESpace> space_b, |
| 19 | // int inda, int indb, |
| 20 | shared_ptr<DifferentialOperator> diffop, shared_ptr<CoefficientFunction> trial_cf, |
| 21 | VorB vb, const Region * reg, LocalHeap & clh, |
| 22 | shared_ptr<BitArray> range_dofs = nullptr, |
| 23 | bool localop = false, bool parmat = true, bool use_simd = true, |
| 24 | int bonus_intorder_ab = 0, int bonus_intorder_bb = 0) |
| 25 | { |
| 26 | static Timer t ("ConvertOperator"); |
| 27 | RegionTimer regt(t); |
| 28 | static Timer tass("CovnertOperator - assemble"); |
| 29 | static Timer telmat("CovnertOperator - elmat"); |
| 30 | |
| 31 | /** Solves gfb = diffop(gfa), where gfb is a function from space_b and gfa one from space_a **/ |
| 32 | |
| 33 | /** Other Sparse-Matrix templates are not instantiated **/ |
| 34 | static_assert( (DIMA <= MAX_SYS_DIM) && (DIMB <= MAX_SYS_DIM) && ( (DIMA == DIMB) || (DIMA == 1) || (DIMB == 1) ) ); |
| 35 | |
| 36 | // /** Compound space is used for element-coloring **/ |
| 37 | // shared_ptr<FESpace> space_a = (*comp_space)[inda]; |
| 38 | // shared_ptr<FESpace> space_b = (*comp_space)[indb]; |
| 39 | auto ma = space_b->GetMeshAccess(); |
| 40 | |
| 41 | if ( parmat && (space_a->IsParallel() != space_b->IsParallel()) ) |
| 42 | { throw Exception("Cannot form ConvertOperator between a parallel and a local space!"); } |
| 43 | |
| 44 | if ( parmat && space_a->IsParallel() && space_b->IsParallel() ) { |
| 45 | NG_MPI_Comm comma = space_a->GetParallelDofs()->GetCommunicator(), commb = space_b->GetParallelDofs()->GetCommunicator(); |
| 46 | if (comma != commb) |
| 47 | { throw Exception("Cannot convert between spaces defined on different Communicators!"); } |
| 48 | } |
| 49 | |
| 50 | /** Proxies and Integrators **/ |
| 51 | |
| 52 | shared_ptr<CoefficientFunction> trial_a; |
| 53 | if (trial_cf != nullptr) |
| 54 | { trial_a = trial_cf; } |
| 55 | else if ( diffop != nullptr ) { |
| 56 | trial_a = make_shared<ProxyFunction>(space_a, false, false, diffop, |
| 57 | nullptr, nullptr, nullptr, nullptr, nullptr); |
| 58 | } |
| 59 | else { // Probably most of the time |
| 60 | trial_a = make_shared<ProxyFunction>(space_a, false, false, space_a->GetEvaluator(vb), |
| 61 | nullptr, nullptr, nullptr, nullptr, nullptr); |
| 62 | } |
| 63 | |
| 64 | auto trial_b = make_shared<ProxyFunction>(space_b, false, false, space_b->GetEvaluator(vb), |
| 65 | nullptr, nullptr, nullptr, nullptr, nullptr); |
| 66 | |
| 67 | if (!space_b->GetAdditionalEvaluators().Used("dual")) |
| 68 | throw Exception(string("Dual diffop does not exist for ") + space_b->GetClassName() + string("!")); |
| 69 | auto dual_evaluator = space_b->GetAdditionalEvaluators()["dual"]; |
| 70 | for (VorB avb = dual_evaluator->VB(); avb < vb; avb++) { |
| 71 | dual_evaluator = dual_evaluator->GetTrace(); |
| 72 | if ( dual_evaluator == nullptr ) |
| 73 | { throw Exception(space_b->GetClassName() + string(" has no dual trace operator for vb = ") |
| 74 | + to_string(avb) + string(" -> ") + to_string(avb + 1) + string("!")); } |
| 75 | } |
no test coverage detected