codim-2 divshape: on edge + face jumps
| 122 | |
| 123 | // codim-2 divshape: on edge + face jumps |
| 124 | virtual void CalcDivShapeJump(const IntegrationPoint &ip, SliceMatrix<double> div_e_shapes) const |
| 125 | { |
| 126 | if (ip.VB() != BBND) // VB() is the enum of the integration point = {VOL, BND, BBND, BBBND} |
| 127 | throw Exception("only BBoundary integration points allowed"); |
| 128 | |
| 129 | if (Dim() == 3) |
| 130 | { |
| 131 | int enr = ip.FacetNr(); |
| 132 | IVec<2> nfaces = e2f_oriented()[enr]; |
| 133 | |
| 134 | div_e_shapes = 0; |
| 135 | |
| 136 | for (int j = 0; j < 2; j++) |
| 137 | { |
| 138 | IntegrationPoint ip_face = ip; |
| 139 | ip_face.SetFacetNr(nfaces[j], BND); |
| 140 | |
| 141 | Matrix<double> shape_f(ndof, 9); |
| 142 | CalcShape(ip_face, shape_f); |
| 143 | |
| 144 | Vec<3> t = GetTangents()[enr]; |
| 145 | Vec<3> n = ElementTopology::GetNormals<3>(ET_TET)[nfaces[j]]; |
| 146 | |
| 147 | Vec<3> nu = Cross(n, t); |
| 148 | |
| 149 | double sign = (j == 0) ? -1 : 1; |
| 150 | |
| 151 | double factor= 1.0/sqr(L2Norm(n)); |
| 152 | factor *= sign; |
| 153 | for (int k = 0; k < ndof; k++) |
| 154 | { |
| 155 | FlatMatrix<> shapei = shape_f.Row(k).AsMatrix(3,3); |
| 156 | Vec<3> sigma_nu = shapei * nu; |
| 157 | div_e_shapes.Row(k) += factor * sigma_nu; |
| 158 | } |
| 159 | } |
| 160 | } |
| 161 | else |
| 162 | { |
| 163 | int vnr = ip.FacetNr(); |
| 164 | IVec<2> nedges = v2e_oriented()[vnr]; |
| 165 | |
| 166 | div_e_shapes = 0; |
| 167 | |
| 168 | for (int j = 0; j < 2; j++) |
| 169 | { |
| 170 | IntegrationPoint ip_face = ip; |
| 171 | ip_face.SetFacetNr(nedges[j], BND); |
| 172 | |
| 173 | Matrix<double> shape_f(ndof, 4); |
| 174 | CalcShape(ip_face, shape_f); |
| 175 | |
| 176 | Vec<2> n = ElementTopology::GetNormals<2>(ET_TRIG)[nedges[j]]; |
| 177 | double sign = (j == 0) ? -1 : 1; |
| 178 | double factor = 1/sqr(L2Norm(n)); |
| 179 | Vec<2> tau(-n(1),n(0)); |
| 180 | |
| 181 | for (int k = 0; k < ndof; k++) |
no test coverage detected