| 798 | } |
| 799 | |
| 800 | virtual void CalcFacetVector (const FiniteElement & volumefel, int LocalFacetNr, |
| 801 | const ElementTransformation & eltrans, FlatArray<int> & ElVertices, |
| 802 | const ElementTransformation & seltrans, |
| 803 | FlatVector<double> elvec, LocalHeap & lh) const |
| 804 | { |
| 805 | static int timer = NgProfiler::CreateTimer ("DGFacet_NeumannBoundaryIntegrator"); |
| 806 | |
| 807 | NgProfiler::RegionTimer reg (timer); |
| 808 | const ScalarFiniteElement<D> * fel1_l2 = |
| 809 | dynamic_cast<const ScalarFiniteElement<D>*> (&volumefel); |
| 810 | ELEMENT_TYPE eltype1 = volumefel.ElementType(); |
| 811 | int nd1 = fel1_l2->GetNDof(); |
| 812 | elvec = 0.0; |
| 813 | |
| 814 | FlatVector<> mat1_shape(nd1, lh); |
| 815 | FlatVector<> mat1_dudn(nd1, lh); |
| 816 | |
| 817 | Facet2ElementTrafo transform1(eltype1);//,ElVertices); at domain boundaries don't change orientation: orientation should still coincide with the orientation of the surface element transformation |
| 818 | |
| 819 | const NORMAL * normals1 = ElementTopology::GetNormals(eltype1); |
| 820 | |
| 821 | HeapReset hr(lh); |
| 822 | ELEMENT_TYPE etfacet = ElementTopology::GetFacetType (eltype1, LocalFacetNr); |
| 823 | |
| 824 | Vec<D> normal_ref1; |
| 825 | for (int i=0; i<D; i++){ |
| 826 | normal_ref1(i) = normals1[LocalFacetNr][i]; |
| 827 | } |
| 828 | int maxorder = fel1_l2->Order(); |
| 829 | const IntegrationRule & ir_facet = |
| 830 | SelectIntegrationRule (etfacet, 2*maxorder); |
| 831 | if (maxorder==0) maxorder=1; |
| 832 | |
| 833 | for (int l = 0; l < ir_facet.GetNIP(); l++) |
| 834 | { |
| 835 | IntegrationPoint ip1 = transform1(LocalFacetNr, ir_facet[l]); |
| 836 | |
| 837 | MappedIntegrationPoint<D,D> sip1 (ip1, eltrans); |
| 838 | double lam = coef_lam->Evaluate(sip1); |
| 839 | |
| 840 | MappedIntegrationPoint<D-1,D> sips (ir_facet[l], seltrans); |
| 841 | |
| 842 | // Mat<D> jac1 = sip1.GetJacobian(); |
| 843 | Mat<D> inv_jac1 = sip1.GetJacobianInverse(); |
| 844 | double det1 = sip1.GetJacobiDet(); |
| 845 | |
| 846 | Vec<D> normal1 = det1 * Trans (inv_jac1) * normal_ref1; |
| 847 | double len1 = L2Norm (normal1); |
| 848 | normal1 /= len1; |
| 849 | |
| 850 | fel1_l2->CalcShape(sip1.IP(), mat1_shape); |
| 851 | |
| 852 | // Vec<D> invjac_normal1 = inv_jac1 * normal1; |
| 853 | |
| 854 | double fac = len1*ir_facet[l].Weight()*lam; |
| 855 | elvec += fac * mat1_shape; |
| 856 | } |
| 857 | } |
no test coverage detected