| 5139 | #ifdef OLD |
| 5140 | template <int DIM> |
| 5141 | void VectorL2FESpace :: |
| 5142 | SolveMPiola (CoefficientFunction * rho, BaseVector & vec, |
| 5143 | LocalHeap & lh) const |
| 5144 | { |
| 5145 | static Timer t("SolveM - Piola"); RegionTimer reg(t); |
| 5146 | |
| 5147 | IterateElements |
| 5148 | (*this, VOL, lh, |
| 5149 | [&rho, &vec,this] (FESpace::Element el, LocalHeap & lh) |
| 5150 | { |
| 5151 | auto & fel = static_cast<const VectorFiniteElement&>(el.GetFE()); |
| 5152 | auto & feli = static_cast<const BaseScalarFiniteElement&>(fel[0]); |
| 5153 | const ElementTransformation & trafo = el.GetTrafo(); |
| 5154 | |
| 5155 | Array<int> dnums(fel.GetNDof(), lh); |
| 5156 | GetDofNrs (el.Nr(), dnums); |
| 5157 | |
| 5158 | FlatVector<double> elx(feli.GetNDof()*DIM, lh); |
| 5159 | vec.GetIndirect(dnums, elx); |
| 5160 | auto melx = elx.AsMatrix(DIM, feli.GetNDof()); |
| 5161 | |
| 5162 | FlatVector<double> diag_mass(feli.GetNDof(), lh); |
| 5163 | feli.GetDiagMassMatrix (diag_mass); |
| 5164 | |
| 5165 | bool curved = trafo.IsCurvedElement(); |
| 5166 | if (rho && !rho->ElementwiseConstant()) curved = true; |
| 5167 | curved = false; // curved not implemented |
| 5168 | |
| 5169 | if (!curved) |
| 5170 | { |
| 5171 | IntegrationRule ir(fel.ElementType(), 0); |
| 5172 | BaseMappedIntegrationRule & mir = trafo(ir, lh); |
| 5173 | Mat<DIM> trans = (1/mir[0].GetMeasure()) * Trans(mir[0].GetJacobian()) * mir[0].GetJacobian(); |
| 5174 | Mat<DIM> invtrans = Inv(trans); |
| 5175 | |
| 5176 | // double jac = mir[0].GetMeasure(); |
| 5177 | // if (rho) jac *= rho->Evaluate(mir[0]); |
| 5178 | // diag_mass *= jac; |
| 5179 | |
| 5180 | for (int i = 0; i < melx.Width(); i++) |
| 5181 | { |
| 5182 | Vec<DIM> hv = melx.Col(i); |
| 5183 | hv /= diag_mass(i); |
| 5184 | melx.Col(i) = invtrans * hv; |
| 5185 | } |
| 5186 | } |
| 5187 | /* |
| 5188 | else |
| 5189 | { |
| 5190 | SIMD_IntegrationRule ir(fel.ElementType(), 2*fel.Order()); |
| 5191 | auto & mir = trafo(ir, lh); |
| 5192 | FlatVector<SIMD<double>> pntvals(ir.Size(), lh); |
| 5193 | FlatMatrix<SIMD<double>> rhovals(1, ir.Size(), lh); |
| 5194 | if (rho) rho->Evaluate (mir, rhovals); |
| 5195 | |
| 5196 | for (int i = 0; i < melx.Height(); i++) |
| 5197 | melx.Row(i) /= diag_mass(i); |
| 5198 | for (int comp = 0; comp < dimension; comp++) |
nothing calls this directly
no test coverage detected