31 bool point_in_convex(
const geometric_trans &geoTrans,
36 return (geoTrans.convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
41 for (
auto &coord : x) {
42 if (coord < 0.0) coord = 0.0;
43 if (coord > 1.0) coord = 1.0;
47 auto pgt_torus = std::dynamic_pointer_cast<const torus_geom_trans>(pgt);
49 orig_pgt = pgt_torus ? pgt_torus->get_original_transformation()
53 auto nb_simplices = pbasic_convex_ref->simplexified_convex()->nb_convex();
55 if (nb_simplices == 1) {
56 auto sum_coordinates = 0.0;
57 for (
const auto &coord : x) sum_coordinates += coord;
59 if (sum_coordinates > 1.0) gmm::scale(x, 1.0 / sum_coordinates);
61 else if (pgt->dim() == 3 && nb_simplices != 4) {
62 auto sum_coordinates = x[0] + x[1];
63 if (sum_coordinates > 1.0) {
64 x[0] /= sum_coordinates;
65 x[1] /= sum_coordinates;
72 bool project_into_element) {
73 bool converged =
true;
74 return invert(n, n_ref, converged, IN_EPS, project_into_element);
80 bool project_into_element) {
82 n_ref.resize(pgt->structure()->dim());
85 return invert_lin(n, n_ref, IN_EPS);
87 return invert_nonlin(n, n_ref, IN_EPS, converged,
false,
88 project_into_element);
95 mult(transposed(B), y, n_ref);
96 y = pgt->transform(n_ref, G);
97 add(gmm::scaled(n, -1.0), y);
99 return point_in_convex(*pgt, n_ref, gmm::vect_norm2(y), IN_EPS);
102 void geotrans_inv_convex::update_B() {
104 pgt->compute_K_matrix(G, pc, K);
105 gmm::mult(gmm::transposed(K), K, CS);
106 bgeot::lu_inverse(&(*(CS.begin())), P);
112 base_matrix KT(K.nrows(), K.ncols());
113 pgt->compute_K_matrix(G, pc, KT);
114 gmm::copy(gmm::transposed(KT), K);
116 bgeot::lu_inverse(&(*(K.begin())), P); B.swap(K);
120 class geotrans_inv_convex_bfgs {
121 geotrans_inv_convex &gic;
124 geotrans_inv_convex_bfgs(geotrans_inv_convex &gic_,
125 const base_node &xr) : gic(gic_), xreal(xr) {}
126 scalar_type operator()(
const base_node& x)
const {
127 base_node r = gic.pgt->transform(x, gic.G) - xreal;
130 void operator()(
const base_node& x, base_small_vector& gr)
const {
131 gic.pgt->poly_vector_grad(x, gic.pc);
133 base_node r = gic.pgt->transform(x, gic.G) - xreal;
135 gmm::mult(gmm::transposed(gic.K), r, gr);
139 void geotrans_inv_convex::update_linearization() {
141 const convex_ind_ct &dir_pt_ind = pgt->structure()->ind_dir_points();
142 const stored_point_tab &ref_nodes = pgt->geometric_nodes();
144 has_linearized_approx =
true;
146 auto n_points = dir_pt_ind.size();
147 auto N_ref = ref_nodes.begin()->size();
149 std::vector<base_node> dir_pts, dir_pts_ref;
150 for (
auto i : dir_pt_ind) {
151 dir_pts.push_back(base_node(N));
152 gmm::copy(mat_col(G, i), dir_pts.back());
153 dir_pts_ref.push_back(ref_nodes[i]);
156 base_matrix K_lin(N, n_points - 1),
157 B_transp_lin(n_points - 1, N),
158 K_ref_lin(N_ref, n_points - 1);
161 P_ref_lin = dir_pts_ref[0];
163 for (
size_type i = 1; i < n_points; ++i) {
164 add(dir_pts[i], gmm::scaled(P_lin, -1.0), mat_col(K_lin, i - 1));
165 add(dir_pts_ref[i], gmm::scaled(P_ref_lin, -1.0),
166 mat_col(K_ref_lin, i - 1));
169 if (K_lin.nrows() == K_lin.ncols()) {
171 gmm::copy(K_lin, B_transp_lin);
174 base_matrix temp(n_points - 1, n_points - 1);
175 mult(transposed(K_lin), K_lin, temp);
177 mult(temp, transposed(K_lin), B_transp_lin);
180 K_ref_B_transp_lin.base_resize(N_ref, N);
181 mult(K_ref_lin, B_transp_lin, K_ref_B_transp_lin);
188 bool geotrans_inv_convex::invert_nonlin(
const base_node& xreal,
189 base_node& x, scalar_type IN_EPS,
192 bool project_into_element) {
195 base_node x0_ref(P), diff(N);
198 x0_ref = pgt->geometric_nodes()[0];
200 for (
size_type j = 1; j < pgt->nb_points(); ++j) {
204 x0_ref = pgt->geometric_nodes()[j];
208 scalar_type res0 = std::numeric_limits<scalar_type>::max();
209 if (has_linearized_approx) {
211 add(xreal, gmm::scaled(P_lin, -1.0), diff);
212 mult(K_ref_B_transp_lin, diff, x);
213 gmm::add(P_ref_lin, x);
215 if (project_into_element) project_into_convex(x, pgt);
219 if (res < res0) gmm::copy(x0_ref, x);
224 add(pgt->transform(x, G), gmm::scaled(xreal, -1.0), diff);
226 auto res0 = std::numeric_limits<scalar_type>::max();
229 base_node x0_real(N);
230 while (res > IN_EPS) {
231 if ((gmm::abs(res - res0) < IN_EPS) || (factor < IN_EPS)) {
233 return point_in_convex(*pgt, x, res, IN_EPS);
236 add(gmm::scaled(x0_ref, factor), x);
237 x0_real = pgt->transform(x, G);
238 add(x0_real, gmm::scaled(xreal, -1.0), diff);
242 if (factor < 1.0-IN_EPS) factor = 2.0;
245 pgt->poly_vector_grad(x, pc);
247 mult(transposed(B), diff, x0_ref);
248 add(gmm::scaled(x0_ref, -1.0 * factor), x);
249 if (project_into_element) project_into_convex(x, pgt);
250 x0_real = pgt->transform(x, G);
251 add(x0_real, gmm::scaled(xreal, -1.0), diff);
255 return point_in_convex(*pgt, x, res, IN_EPS);