| 4791 | |
| 4792 | template <int DIM> |
| 4793 | void VectorL2FESpace :: |
| 4794 | SolveM_Dim (CoefficientFunction * rho, BaseVector & vec, Region * def, |
| 4795 | LocalHeap & lh) const |
| 4796 | { |
| 4797 | static Timer t("SolveM - Vec"); RegionTimer reg(t); |
| 4798 | IterateElements |
| 4799 | (*this, VOL, lh, |
| 4800 | [&rho, &vec, def, this] (FESpace::Element el, LocalHeap & lh) |
| 4801 | { |
| 4802 | auto & fel = static_cast<const VectorFiniteElement&>(el.GetFE()); |
| 4803 | auto & feli = static_cast<const BaseScalarFiniteElement&>(fel[0]); |
| 4804 | const ElementTransformation & trafo = el.GetTrafo(); |
| 4805 | |
| 4806 | Array<int> dnums(fel.GetNDof(), lh); |
| 4807 | GetDofNrs (el.Nr(), dnums); |
| 4808 | |
| 4809 | FlatVector<double> elx(feli.GetNDof()*DIM, lh); |
| 4810 | |
| 4811 | if (def && !def->Mask()[ma->GetElIndex(el)]) |
| 4812 | { |
| 4813 | elx = 0.0; |
| 4814 | vec.SetIndirect (dnums, elx); |
| 4815 | return; |
| 4816 | } |
| 4817 | |
| 4818 | vec.GetIndirect(dnums, elx); |
| 4819 | auto melx = elx.AsMatrix(DIM, feli.GetNDof()); |
| 4820 | |
| 4821 | FlatVector<double> diag_mass(feli.GetNDof(), lh); |
| 4822 | feli.GetDiagMassMatrix (diag_mass); |
| 4823 | |
| 4824 | bool curved = trafo.IsCurvedElement(); |
| 4825 | if (rho && !rho->ElementwiseConstant()) curved = true; |
| 4826 | // curved = false; // curved not implemented |
| 4827 | |
| 4828 | if (!curved) |
| 4829 | { |
| 4830 | IntegrationRule ir(fel.ElementType(), 0); |
| 4831 | MappedIntegrationRule<DIM,DIM> mir(ir, trafo, lh); |
| 4832 | |
| 4833 | Mat<DIM,DIM> rhoi(0.0); |
| 4834 | if (!rho) |
| 4835 | rhoi = Identity(DIM); |
| 4836 | else if (rho->Dimension() == 1) |
| 4837 | rhoi = rho->Evaluate(mir[0]) * Identity(DIM); |
| 4838 | else |
| 4839 | rho -> Evaluate(mir[0], FlatVector<> (DIM*DIM, &rhoi(0,0))); |
| 4840 | |
| 4841 | Mat<DIM> trans(0.0); |
| 4842 | if (piola) |
| 4843 | trans = (1/mir[0].GetMeasure()) * Trans(mir[0].GetJacobian()) * rhoi * mir[0].GetJacobian(); |
| 4844 | else if (covariant) |
| 4845 | trans = mir[0].GetMeasure() * mir[0].GetJacobianInverse() * rhoi * Trans(mir[0].GetJacobianInverse()); |
| 4846 | else |
| 4847 | trans = mir[0].GetMeasure() * rhoi; |
| 4848 | |
| 4849 | Mat<DIM> invtrans = Inv(trans); |
| 4850 |
nothing calls this directly
no test coverage detected