FastMarchingInterface.h 6.92 KB
Newer Older
1 2 3 4 5
/* 
    This file is a part of eXlibris C++ Library
    under the GNU General Public License:
    See the LICENSE.md files for terms and 
    conditions.
6
*/
7

Kevin Moreau's avatar
Kevin Moreau committed
8 9
#ifndef _FastMarchingInterface_h_
#define _FastMarchingInterface_h_
10 11


Kevin Moreau's avatar
Kevin Moreau committed
12 13 14 15 16
// Xfem
#include "xRegion.h"
//
#ifdef HAVE_FASTMARCHING
#include "meshinterfacexRegion.h"
Kévin Moreau's avatar
Kévin Moreau committed
17
#include "FM.h"
18 19 20
#include "FMEntityStorage.h"
template < typename T >
using datamanagerFM_t = xfastmarching::entitystorage < xfastmarching::meshinterfacexRegion, AOMD::mVertex,T >;
Kevin Moreau's avatar
Kevin Moreau committed
21 22 23 24
#endif
// std
#include <functional>

Kevin Moreau's avatar
Kevin Moreau committed
25 26 27 28 29 30 31 32 33 34
class modesValues{
public:
  modesValues(){};
  modesValues(size_t n, double val):values(n,val){};
  modesValues(const modesValues &in):values(in.values){};
  double & operator[](const size_t &i){ return values[i];}
  modesValues & operator*=(const double &scal){
    for_each(values.begin(), values.end(), [&scal](double & val){ val*= scal;} );
    return *this;
  }
35 36 37 38 39
  modesValues & operator/=(const double &scal){
    for_each(values.begin(), values.end(), [&scal](double & val){ val/= scal;} );
    return *this;
  }
  
Kevin Moreau's avatar
Kevin Moreau committed
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
  modesValues & operator+=(const modesValues &in){
    transform( values.begin(), values.end(), in.values.begin(), values.begin(), [] (double val1, double val2 ){return val1+val2;});
    return *this;
  }
  modesValues& operator=(const modesValues &in) {
    values=in.values;
    return *this;
  }
  friend modesValues operator+(modesValues lhs,  const modesValues& rhs)
  {
    return lhs += rhs; // reuse compound assignment and return the result by value
  }
  friend modesValues operator*(const double &lhs, modesValues rhs){
    return rhs *= lhs;
  }
  friend modesValues operator*(modesValues rhs, const double &lhs){
    return rhs *= lhs;
  }
58 59 60 61
  friend modesValues operator/(modesValues rhs, const double &lhs){
    return rhs /= lhs;
  }
  
Kevin Moreau's avatar
Kevin Moreau committed
62 63 64 65
  std::vector<double> values;
private :
};

66 67 68 69
template< typename ITERKNOWN
         ,typename ITERKNOWNV
         ,typename ITERTOFIND
         ,typename ITERKNOWNT>
Kévin Moreau's avatar
Kévin Moreau committed
70 71
void FastMarchingReinit(std::function<double (AOMD::mVertex*)> get_val,
                        std::function<void (AOMD::mVertex*, double)> set_val,
Kevin Moreau's avatar
Kevin Moreau committed
72
                        const xfem::xRegion& region,
73 74
                        ITERKNOWN known_begin, ITERKNOWN known_end,
                        ITERKNOWNV known_volatil_begin, ITERKNOWNV known_volatil_end,
Kévin Moreau's avatar
Kévin Moreau committed
75
                        ITERTOFIND tofind_begin, ITERTOFIND tofind_end,
76
                        ITERKNOWNT known_it,
77
                        double epsilon_ratio,
Kévin Moreau's avatar
Kévin Moreau committed
78
                        std::function<double (const AOMD::mVertex&)> f_func=[](const AOMD::mVertex&){return 1.;}) {
Kevin Moreau's avatar
Kevin Moreau committed
79 80 81 82
#ifndef HAVE_FASTMARCHING
  std::cout << "FastMarchingReinit requires HAVE_FASTMARCHING defined to work" << std::endl;
  throw;
#else
83 84 85
  xfastmarching::meshinterfacexRegion mi(region);
  xfastmarching::entitystorage<xfastmarching::meshinterfacexRegion, AOMD::mVertex, double> ls(mi);
  xfastmarching::entitystorage<xfastmarching::meshinterfacexRegion, AOMD::mVertex, xfastmarching::vector3d<double> > gls(mi);
86 87 88 89 90
  for(auto it=known_begin; it!=known_end; ++it) {
    auto v=static_cast<AOMD::mVertex*>(*it);
    ls.set(*v, get_val(v));
  }
  for(auto it=known_volatil_begin; it!=known_volatil_end; ++it) {
Kévin Moreau's avatar
Kévin Moreau committed
91 92
    auto v=static_cast<AOMD::mVertex*>(*it);
    ls.set(*v, get_val(v));
Kevin Moreau's avatar
Kevin Moreau committed
93
  }
Kévin Moreau's avatar
Kévin Moreau committed
94 95 96
  for(auto it=tofind_begin; it!=tofind_end; ++it) {
    ls.set(*static_cast<AOMD::mVertex*>(*it), std::numeric_limits<double>::max());
  }
97 98 99 100 101
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWN> trial_known_conv_begin(known_begin);
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWN> trial_known_conv_end(known_end);
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWNV> trial_known_volatil_conv_begin(known_volatil_begin);
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWNV> trial_known_volatil_conv_end(known_volatil_end);
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWNT> known_conv_it(known_it);
102
  xfastmarching::fmeik<datamanagerFM_t>(mi, ls, trial_known_conv_begin, trial_known_conv_end,trial_known_volatil_conv_begin, trial_known_volatil_conv_end, known_conv_it, f_func, epsilon_ratio, gls);
Kévin Moreau's avatar
Kévin Moreau committed
103 104
  for(auto it=region.begin(0); it!=region.end(0); ++it) {
    auto v=static_cast<AOMD::mVertex*>(*it);
Kevin Moreau's avatar
Kevin Moreau committed
105
    double val;
Kévin Moreau's avatar
Kévin Moreau committed
106 107
    ls.get(*v, val);
    set_val(v, val);
Kevin Moreau's avatar
Kevin Moreau committed
108 109 110 111
  }
#endif
}

Kévin Moreau's avatar
Kévin Moreau committed
112
template<typename ITERKNOWN,
113 114
         typename ITERTOFIND,
         typename MATRIX>
Kévin Moreau's avatar
Kévin Moreau committed
115 116
void FastMarchingModeExtension(std::function<double (AOMD::mVertex*)> get_val,
                               std::function<void (AOMD::mVertex*, double)> set_val,
Kevin Moreau's avatar
Kevin Moreau committed
117
                               const xfem::xRegion& region,
Kévin Moreau's avatar
Kévin Moreau committed
118 119
                               ITERKNOWN known_begin, ITERKNOWN known_end,
                               ITERTOFIND tofind_begin, ITERTOFIND tofind_end,
120 121
                               MATRIX& coeffs,
                               double epsilon_ratio,
Kévin Moreau's avatar
Kévin Moreau committed
122
                               const std::function<double (const AOMD::mVertex&)>& f_func=[](const AOMD::mVertex&){return 1.;}) {
Kevin Moreau's avatar
Kevin Moreau committed
123 124 125 126
#ifndef HAVE_FASTMARCHING
  std::cout << "FastMarchingModeExtension requires HAVE_FASTMARCHING defined to work" << std::endl;
  throw;
#else
127 128 129 130
  xfastmarching::meshinterfacexRegion mi(region);
  xfastmarching::entitystorage<xfastmarching::meshinterfacexRegion, AOMD::mVertex, double> ls(mi);
  xfastmarching::entitystorage<xfastmarching::meshinterfacexRegion, AOMD::mVertex, xfastmarching::vector3d<double> > gls(mi);
  xfastmarching::entitystorage<xfastmarching::meshinterfacexRegion, AOMD::mVertex, modesValues> vn(mi);
Kévin Moreau's avatar
Kévin Moreau committed
131
  const int nb_modes=coeffs.nline();
Kevin Moreau's avatar
Kevin Moreau committed
132
  int i=0;
Kévin Moreau's avatar
Kévin Moreau committed
133 134 135
  for(auto it=known_begin; it!=known_end; ++it, ++i) {
    auto v=static_cast<AOMD::mVertex*>(*it);
    ls.set(*v, get_val(v));
Kevin Moreau's avatar
Kevin Moreau committed
136 137
    modesValues mv(nb_modes, 0.);
    mv[i]=1.;
Kévin Moreau's avatar
Kévin Moreau committed
138 139 140 141
    vn.set(*v, mv);
  }
  for(auto it=tofind_begin; it!=tofind_end; ++it) {
    ls.set(*static_cast<AOMD::mVertex*>(*it), std::numeric_limits<double>::max());
Kevin Moreau's avatar
Kevin Moreau committed
142
  }
143
  const std::set<const AOMD::mVertex*> empty;
144 145 146 147
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWN> known_conv_begin(known_begin);
  xfastmarching::entitytovertexiteratorconvertor<ITERKNOWN> known_conv_end(known_end);
  xfastmarching::entitytovertexiteratorconvertor< std::set<const AOMD::mVertex*>::const_iterator> known_volatil_conv_begin(empty.begin());
  xfastmarching::entitytovertexiteratorconvertor< std::set<const AOMD::mVertex*>::const_iterator> known_volatil_conv_end(empty.end());
Kévin Moreau's avatar
Kévin Moreau committed
148 149
  std::vector<const AOMD::mVertex*> vec;
  vec.reserve(nb_modes);
150
  xfastmarching::entitytovertexiteratorconvertor<std::back_insert_iterator<std::vector<const AOMD::mVertex*>>> dummy_it(std::back_inserter(vec));
151
  xfastmarching::fmeik<datamanagerFM_t>(mi, ls, known_conv_begin, known_conv_end,known_volatil_conv_begin, known_volatil_conv_end, dummy_it, f_func, epsilon_ratio, gls, vn);
Kevin Moreau's avatar
Kevin Moreau committed
152 153
  int j=0;
  for(auto it=region.begin(0); it!=region.end(0); ++it, ++j) {
Kévin Moreau's avatar
Kévin Moreau committed
154 155 156 157
    auto v=static_cast<AOMD::mVertex*>(*it);
    double val;
    ls.get(*v, val);
    set_val(v, val);
Kevin Moreau's avatar
Kevin Moreau committed
158
    modesValues mv;
Kévin Moreau's avatar
Kévin Moreau committed
159
    vn.get(*v, mv);
Kevin Moreau's avatar
Kevin Moreau committed
160 161 162 163 164 165
    for(int i=0; i<nb_modes; ++i) {
      coeffs(i,j)=mv[i];
    }
  }
#endif
}
Kévin Moreau's avatar
Kévin Moreau committed
166

Kevin Moreau's avatar
Kevin Moreau committed
167
#endif