| 4932 | |
| 4933 | template <int DIM> |
| 4934 | shared_ptr<BaseMatrix> |
| 4935 | VectorL2FESpace :: CreateMassOperator_Dim (shared_ptr<CoefficientFunction> rho, |
| 4936 | shared_ptr<Region> defon, |
| 4937 | bool inverse, |
| 4938 | LocalHeap & lh) const |
| 4939 | { |
| 4940 | shared_ptr<BaseMatrix> sum; |
| 4941 | |
| 4942 | bool optimize_constant = rho->ElementwiseConstant(); |
| 4943 | if (optimize_constant) |
| 4944 | { |
| 4945 | shared_ptr<BaseMatrix> mat; |
| 4946 | |
| 4947 | auto & vfe = static_cast<const VectorFiniteElement&>(GetFE(ElementId(VOL, 0), lh)); |
| 4948 | auto & fe = static_cast<const BaseScalarFiniteElement&>(vfe[0]); |
| 4949 | |
| 4950 | Vector<> diag_mass(fe.GetNDof()); |
| 4951 | dynamic_cast<const BaseScalarFiniteElement&>(fe).GetDiagMassMatrix(diag_mass); |
| 4952 | auto ma = GetMeshAccess(); |
| 4953 | Array<size_t> elindex(ma->GetNE()); |
| 4954 | elindex = size_t(-1); |
| 4955 | size_t count = 0; |
| 4956 | for (const auto & el : Elements(VOL)) |
| 4957 | if(DefinedOn(el)) |
| 4958 | elindex[el.Nr()] = count++; |
| 4959 | Vector<Mat<DIM,DIM>> elscale(count); |
| 4960 | |
| 4961 | IterateElements |
| 4962 | (*this, VOL, lh, |
| 4963 | [this, inverse, &ma, &elscale, &defon,&rho, &elindex] (FESpace::Element el, LocalHeap & lh) |
| 4964 | { |
| 4965 | auto & fel = static_cast<const BaseScalarFiniteElement&>(el.GetFE()); |
| 4966 | const ElementTransformation & trafo = el.GetTrafo(); |
| 4967 | |
| 4968 | IntegrationRule ir(fel.ElementType(), 0); |
| 4969 | BaseMappedIntegrationRule & mir = trafo(ir, lh); |
| 4970 | |
| 4971 | Mat<DIM> transrho; |
| 4972 | if ( (defon && !defon->Mask()[ma->GetElIndex(el)]) || el.is_curved) |
| 4973 | transrho = 0; |
| 4974 | else |
| 4975 | { |
| 4976 | Mat<DIM> rhoi(0.0); |
| 4977 | if (!rho) |
| 4978 | rhoi = Identity(DIM); |
| 4979 | else if (rho->Dimension() == 1) |
| 4980 | rhoi = rho->Evaluate(mir[0]) * Identity(DIM); |
| 4981 | else |
| 4982 | rho -> Evaluate(mir[0], FlatVector<> (DIM*DIM, &rhoi(0,0))); |
| 4983 | |
| 4984 | if (piola) |
| 4985 | transrho = 1/mir[0].GetMeasure() * Trans(mir[0].GetJacobian()) * rhoi * mir[0].GetJacobian(); |
| 4986 | else if (covariant) |
| 4987 | transrho = mir[0].GetMeasure() * Inverse(mir[0].GetJacobian()) * rhoi * Trans(Inverse(mir[0].GetJacobian())); |
| 4988 | else |
| 4989 | transrho = mir[0].GetMeasure() * rhoi; |
| 4990 | |
| 4991 | if (inverse) |
nothing calls this directly
no test coverage detected