| 2344 | |
| 2345 | template <typename FEL, typename MIP, typename MAT> |
| 2346 | static void GenerateMatrix (const FEL & bfel, const MIP & mip, |
| 2347 | MAT && bmat, LocalHeap & lh) |
| 2348 | { |
| 2349 | auto & fel = static_cast<const VectorFiniteElement&> (bfel); |
| 2350 | auto & fel_u = static_cast<const ScalarFiniteElement<DIM_SPACE>&> (fel[0]); |
| 2351 | |
| 2352 | int order = fel_u.Order(); |
| 2353 | auto & trafo = mip.GetTransformation(); |
| 2354 | |
| 2355 | HDivHighOrderFE<ET_TRIG> fel_hdiv(order-1); |
| 2356 | L2HighOrderFE<ET_TRIG> fel_l2(order-2); |
| 2357 | L2HighOrderFE<ET_TRIG> fel_koszul( max(order-3, -1) ); |
| 2358 | // FacetFE<ET_TRIG> fel_facet; |
| 2359 | FacetFE<ET_TRIG> & fel_facet = *new (lh) FacetFE<ET_TRIG>; |
| 2360 | fel_facet.SetOrder(order-1); |
| 2361 | Array<int> vnums = { 1, 2, 3 } ; |
| 2362 | fel_facet.SetVertexNumbers(vnums); |
| 2363 | fel_facet.ComputeNDof(); |
| 2364 | |
| 2365 | Matrix<> mat(fel_hdiv.GetNDof()); |
| 2366 | Matrix<> rhs(fel_hdiv.GetNDof(), DIM_SPACE*fel_u.GetNDof()); |
| 2367 | |
| 2368 | auto eltype = trafo.GetElementType(); |
| 2369 | int nfacet = ElementTopology::GetNFacets(eltype); |
| 2370 | Facet2ElementTrafo transform(eltype); |
| 2371 | |
| 2372 | IntRange r_facet(0, fel_facet.GetNDof()); |
| 2373 | IntRange r_div(r_facet.Next(), r_facet.Next() + fel_l2.GetNDof() - 1); |
| 2374 | IntRange r_koszul(r_div.Next(), fel_hdiv.GetNDof()); |
| 2375 | |
| 2376 | /* |
| 2377 | *testout << "r_facet = " << r_facet << endl; |
| 2378 | *testout << "r_div = " << r_div << endl; |
| 2379 | *testout << "r_koszul = " << r_koszul << endl; |
| 2380 | */ |
| 2381 | |
| 2382 | mat = 0; |
| 2383 | rhs = 0; |
| 2384 | |
| 2385 | FlatMatrix<> shape_hdiv(fel_hdiv.GetNDof(), 2, lh); |
| 2386 | FlatVector<> shape_hdivn(fel_hdiv.GetNDof(), lh); |
| 2387 | FlatVector<> shape_u(fel_u.GetNDof(),lh); |
| 2388 | |
| 2389 | // edge moments |
| 2390 | for (int k = 0; k < nfacet; k++) |
| 2391 | { |
| 2392 | IntRange r_facet_k = fel_facet.GetFacetDofs(k); |
| 2393 | FlatVector<> shape_facet(r_facet_k.Size(), lh); |
| 2394 | |
| 2395 | HeapReset hr(lh); |
| 2396 | ngfem::ELEMENT_TYPE etfacet = ElementTopology::GetFacetType (eltype, k); |
| 2397 | |
| 2398 | IntegrationRule ir_facet(etfacet, 2*fel.Order()); |
| 2399 | IntegrationRule & ir_facet_vol = transform(k, ir_facet, lh); |
| 2400 | |
| 2401 | MappedIntegrationRule<DIM_ELEMENT,DIM_SPACE> mir(ir_facet_vol, trafo, lh); |
| 2402 | mir.ComputeNormalsAndMeasure (eltype, k); |
| 2403 |
nothing calls this directly
no test coverage detected