| 1460 | |
| 1461 | |
| 1462 | IntegrationPoint ProjectPointToTriangleReference(Vec<3> x, const ElementTransformation & trafo) |
| 1463 | { |
| 1464 | IntegrationPoint ip(1./3, 1./3); |
| 1465 | for (int j = 0; j < 5; j++) // SQP steps |
| 1466 | { |
| 1467 | MappedIntegrationPoint<2,3> mip(ip, trafo); |
| 1468 | // dist = || x - (mip+Jac*(uv-ip) + 1/2*Hesse(uv-ip, uv-ip)) || |
| 1469 | Mat<3,2> jac = mip.GetJacobian(); |
| 1470 | Vec<2> ipvec { ip(0), ip(1) }; |
| 1471 | auto Hesse = mip.CalcHesse(); |
| 1472 | Vec<3,Vec<2>> Hesseip |
| 1473 | { |
| 1474 | Hesse[0]*ipvec, |
| 1475 | Hesse[1]*ipvec, |
| 1476 | Hesse[2]*ipvec |
| 1477 | }; |
| 1478 | Vec<3> Hesseipip |
| 1479 | { |
| 1480 | InnerProduct(Hesseip(0), ipvec), |
| 1481 | InnerProduct(Hesseip(1), ipvec), |
| 1482 | InnerProduct(Hesseip(2), ipvec) |
| 1483 | }; |
| 1484 | Mat<3,2> jacphip = jac; |
| 1485 | jacphip.Row(0) -= Hesseip(0); |
| 1486 | jacphip.Row(1) -= Hesseip(1); |
| 1487 | jacphip.Row(2) -= Hesseip(2); |
| 1488 | Mat<2,2> a = Trans(jac)*jac; |
| 1489 | Vec<2> b = -Trans(jacphip) * (x-mip.GetPoint()+jac*ipvec + 0.5*Hesseipip); |
| 1490 | Vec<2> uv = MinimizeOnTrig(a, b, 0); |
| 1491 | ip = IntegrationPoint(uv(0), uv(1)); |
| 1492 | } |
| 1493 | return ip; |
| 1494 | } |
| 1495 | |
| 1496 | |
| 1497 | IntegrationRule GetIntegrationRule(Vec<3> x, const ElementTransformation & trafo, int intorder, bool nearfield) |
no test coverage detected