TLSSolverMode.cc 28.9 KB
Newer Older
1
/*
2 3
    This source code is subject to non-permissive licence,
    see the DamageBandDyn/LICENSE file for conditions.
4
*/
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
#include "refactoring.h"
// DamageBandDyn
#include "Algorithm.h"
#include "CFL.h"
#include "Export.h"
#include "Eval.h"
#include "Geom.h"
#include "Observer.h"
#include "SearchForInitiation.h"
#include "TLSSolverMode.h"
// xTLS
#include "xFrontModeExtensionAlgorithm.h"
#ifdef MODEBND
#include "reInitLevelSetFromFrontAndBoundary.h"
#endif
// SolverBase
#include "xLinearSystemSolverSuperLU.h"

using namespace xfem;

25 26 27
void ClearAttachedEntity(xIter it,
                         xIter end,
                         const unsigned int tag)
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
{
  for(; it!=end; ++it) {
    (*it)->deleteData(tag);
  }
}

void SetTagVertices(boost::filter_iterator<xEntityFilter, xIter> it,
                    boost::filter_iterator<xEntityFilter, xIter> end,
                    const unsigned int tag)
{
  for(; it!=end; ++it) {
    AOMD::mEntity* elt = *it;
    for(int i=0; i<elt->size(0); ++i) {
      elt->get(0,i)->attachInt(tag, 1);
    }
  }
}

void ClearTagVertices(boost::filter_iterator<xEntityFilter, xIter> it,
                      boost::filter_iterator<xEntityFilter, xIter> end,
                      const unsigned int tag)
{
  for(; it!=end; ++it) {
    AOMD::mEntity* elt = *it;
    for(int i=0; i<elt->size(0); ++i) {
      elt->get(0,i)->deleteData(tag);
    }
  }
}

Kevin Moreau's avatar
Kevin Moreau committed
58 59 60 61
void BuildBand(boost::filter_iterator<xEntityFilter, xIter> domain_neg_support_it,
               boost::filter_iterator<xEntityFilter, xIter> domain_neg_support_end,
               const double lower, const double upper,
               const xLevelSet& ls,
62 63
               xSubMesh& sub_mesh,
               const unsigned int tag)
64 65
{
  double ls_val;
Kevin Moreau's avatar
Kevin Moreau committed
66 67
  bool is_inside_band;
  bool has_one_vertex_inside_band;
68
  for(; domain_neg_support_it!=domain_neg_support_end; ++domain_neg_support_it) {
Kevin Moreau's avatar
Kevin Moreau committed
69 70
    is_inside_band = true;
    has_one_vertex_inside_band = false;
71 72 73 74
    AOMD::mEntity* elt = *domain_neg_support_it;
    for(int i=0; i<elt->size(0); ++i) {
      AOMD::mVertex* v = static_cast<AOMD::mVertex*>(elt->get(0, i));
      ls_val = ls(v);
Kevin Moreau's avatar
Kevin Moreau committed
75 76 77
      bool buf = lower<=ls_val && ls_val<=upper;
      is_inside_band = is_inside_band && buf;
      has_one_vertex_inside_band = has_one_vertex_inside_band || buf;
78
    }
Kevin Moreau's avatar
Kevin Moreau committed
79
    if(is_inside_band) {
80
      sub_mesh.add(elt);
81
      for(int i=0; i<elt->size(0); ++i) {
82 83 84
        AOMD::mEntity* v = elt->get(0, i);
        sub_mesh.add(v);
        v->attachEntity(tag, elt);
85 86
      }
    }
Kevin Moreau's avatar
Kevin Moreau committed
87
    else if(has_one_vertex_inside_band) {
88
      sub_mesh.add(elt);
89 90 91 92 93 94
    }
  }
}

class PartialShiftLevelSetModifier : public xLevelSetModifier {
 public:
95
  PartialShiftLevelSetModifier(const xEval<double>& eval_, const xSubMesh& sub_mesh_, unsigned int attach_tag_) :
96
    eval(eval_),
97 98
    sub_mesh(sub_mesh_),
    attach_tag(attach_tag_)
99 100 101 102
  {}
  void visit(xLevelSet& ls, xRegion target)
  {
    double shift;
103 104 105 106
    xIter it=sub_mesh.begin(0);
    for(; it!=sub_mesh.end(0); ++it) {
      AOMD::mVertex* v = static_cast<AOMD::mVertex*>(*it);
      AOMD::mEntity* elt = v->getAttachedEntity(attach_tag); 
107 108 109 110 111 112 113 114 115 116
      xGeomElem geom(elt);
      geom.setUVWForXYZ(v->point());
      eval(&geom, &geom, shift);
      if(shift>0.) {
        ls(v)+=shift;
      }
    }
  }
 private:
  const xEval<double>& eval;
117 118
  const xSubMesh& sub_mesh;
  const unsigned int attach_tag;
119 120
};

121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
#ifdef KALTHOFF_HARD_CODED_TEST
class xAcceptOnBox
{
public:
  xAcceptOnBox(const double xmin_, const double xmax_, const double ymin_, const double ymax_) :
    xmin(xmin_), xmax(xmax_), ymin(ymin_), ymax(ymax_) {}
  bool operator()(AOMD::mEntity* e) const
  {
    if (e->getType()==AOMD::mEntity::VERTEX) {
      Trellis_Util::mPoint p = static_cast<AOMD::mVertex*>(e)->point();
      return (xmin < p(0)) && (p(0) < xmax) && (ymin < p(1)) && (p(1) < ymax);
    }
    else if (e->getType()==AOMD::mEntity::EDGE || e->getType()==AOMD::mEntity::TRI) {
      bool val = true;
      for (int i=0; i<e->size(0); ++i) {
	Trellis_Util::mPoint p = static_cast<AOMD::mVertex*>(e->get(0,i))->point();
	val = val && (xmin < p(0)) && (p(0) < xmax) && (ymin < p(1)) && (p(1) < ymax);
      }
      return val;
    }
    return false;
  }
private:
  const double xmin;
  const double xmax;
  const double ymin;
  const double ymax;
};
#endif

151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
TLSSolverMode::TLSSolverMode(Geom& geom_,
                             xLevelSet* ls_geo_, xLevelSet& ls_in_,
                             xLevelSetModifier& ls_update_, xLevelSetModifier& ls_reinit_,
                             SearchForInitiation& search_for_initiation_,
                             const int dim_, const int integration_order_, const int damage_shape_,
                             const double characteristic_length_, const double diffusive_coefficient_,
                             const double supg_stabilization_factor_, const double mode_min_support_length_factor_,
                             xSpaceFactoryBase& space_factory_,
                             const bool compute_dissipation_,
                             Observer& observer_,
                             PostProcessing& post_processing_) :
  geom(geom_),
  ls_geo(ls_geo_), ls_in(ls_in_),
  ls_update(ls_update_), ls_reinit(ls_reinit_),
  search_for_initiation(search_for_initiation_),
  dim(dim_), integration_order(integration_order_),
  lc(characteristic_length_),
  supg_stabilization_factor(supg_stabilization_factor_), mode_min_support_length_factor(mode_min_support_length_factor_),
  min_element_size(0.), max_element_size(0.),
  eval_normalized_grad_ls(ls_in),
  eval_minus_normalized_grad_ls(eval_normalized_grad_ls),
#ifdef VECTORLS
  eval_adim_ls(ls_in, characteristic_length_, boost::bind(geom.getPartionFunction(), _1, _2, geom.getFilter(Geom::DAMAGEBAND)), geom.getFrontMeshOut(), geom.crackExist(), geom.getCrackUpperCreatorFunction()),
#else
  eval_adim_ls(ls_in, characteristic_length_),
#endif
  eval_damage(&eval_adim_ls, (xtls::damageShapeFunctions::shape_t) damage_shape_, geom.getFilter(Geom::DAMAGEBAND), geom.getFilter(Geom::SANEMATTER)),
  eval_damage_derivative(&eval_adim_ls, (xtls::damageShapeFunctions::shape_t) damage_shape_, geom.getFilter(Geom::DAMAGEBAND), geom.getFilter(Geom::SANEMATTER)),
  front_mode_extension(&space_factory_),
  eval_front_velocity(&space_factory_, compute_dissipation_),
  eval_front_velocity_extended(&space_factory_),
  eval_damage_rate(eval_damage_derivative, eval_front_velocity, geom.getFilter(Geom::DAMAGEBAND)),
  cumulate_nonlocal_dissipation(0.),
  to_move(false),
  observer(observer_),
Kevin Moreau's avatar
Kevin Moreau committed
186
  post_processing(post_processing_),
187 188 189
  narrow_band(NULL),
  narrow_band_attach_tag(AOMD::AOMD_Util::Instance()->lookupMeshDataId("narrow_band_attach_entity")),
  damage_band_attach_tag(AOMD::AOMD_Util::Instance()->lookupMeshDataId("damage_band_attach_entity"))
190 191 192 193 194 195 196 197
{
#ifdef KALTHOFF_HARD_CODED_TEST
  xAcceptOnBox accept_on_box_1(0.05, 0.1, 0.0281, 0.1);
  xAcceptOnBox accept_on_box_2(0.0571, 0.1, 0., 0.0258);
  sensors.registerBox(1, accept_on_box_1);
  sensors.registerBox(2, accept_on_box_2);
#endif
}
198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268

void TLSSolverMode::initialize(const IterDomain& begin_, const IterDomain& end_, const xIntegrationRule& integration_rule_)
{
  std::pair<double, double> min_max_element_size(std::numeric_limits<double>::max(), std::numeric_limits<double>::min());
  ComputeMinMaxInsideRadius compute_min_max_inside_radius;
  Accumulate(compute_min_max_inside_radius, begin_, end_, min_max_element_size);
  min_element_size = min_max_element_size.first;
  max_element_size = min_max_element_size.second;
}

void TLSSolverMode::restore()
{
  std::cout << red("STATE RESTORATION in TLSSolver") << std::endl;
  RestoreLevelSet("ls_in.restore", ls_in);
  RestoreValue("cumulate_nonlocal_dissipation.restore", cumulate_nonlocal_dissipation);
  search_for_initiation.restore();
  updateDomain();
}

void TLSSolverMode::update(const xEval<double>& eval_energy_release_rate_, const xEval<double>& eval_critical_energy_release_rate_,
                           EvalWrapper eval_damage_rate_wrapper_,
                           double time_, int step_, double time_step_,
                           const IterDomain& begin_, const IterDomain& end_, const xIntegrationRule& integration_rule_)
{
  post_processing.exportOnSpace("energy_release_rate", step_, eval_energy_release_rate_, integration_rule_, begin_, end_);

  std::map<Trellis_Util::mPoint, double> hot_points;
  if (search_for_initiation.doSearch())
  {
    xEvalConstant<double> eval_radius(4*max_element_size);
    LocateHotPointsCommand locate_hot_points(eval_energy_release_rate_, eval_critical_energy_release_rate_, eval_radius, hot_points, 0.);
    const IterDomain& search_begin = geom.begin(Geom::SEARCH, dim);
    const IterDomain& search_end = geom.end(Geom::SEARCH, dim);
    const xIntegrationRule& search_integration_rule = *geom.getIntegrationRule(Geom::SEARCH, 0);
    ApplyCommandOnIntegrationRule(locate_hot_points, search_integration_rule, search_begin, search_end);
  }
  search_for_initiation.compute(hot_points);
  ls_update.visit(ls_in, ls_in.getSupport());
}

void TLSSolverMode::compute(xEval<double>& eval_energy_release_rate_, const xEval<double>& eval_critical_energy_release_rate_,
                            EvalWrapper eval_damage_rate_wrapper_,
                            double time_, int step_, double time_step_,
                            const IterDomain& begin_, const IterDomain& end_, const xIntegrationRule& integration_rule_)
{
  if (search_for_initiation.doUpdate()) updateDomain();
  if (search_for_initiation.isFirstInitiated())
  {
    const IterDomain& damage_band_begin = geom.begin(Geom::DAMAGEBANDSUPPORT, dim);
    const IterDomain& damage_band_end   = geom.end  (Geom::DAMAGEBANDSUPPORT, dim);
    const IterDomain& begin             = geom.begin();
    const IterDomain& end               = geom.end();
    const xIntegrationRule& integration_rule = geom.getIntegrationRule();

    post_processing.exportOnSpace("damage_derivative", step_, eval_damage_derivative, integration_rule, begin, end);

    // when predamaged, it avoids to compute the SUPG if not needed
    bool do_supg_computation = true;
    if (geom.getFrontAdvance()==0) {
      std::map<Trellis_Util::mPoint, double> hot_points;
      xEvalConstant<double> eval_radius(4*max_element_size);
      LocateHotPointsCommand locate_hot_points(eval_energy_release_rate_, eval_critical_energy_release_rate_, eval_radius, hot_points, 0.);
      const xIntegrationRule& damage_band_integration_rule = *geom.getIntegrationRule(Geom::DAMAGEBAND, 0);
      ApplyCommandOnIntegrationRule(locate_hot_points, damage_band_integration_rule, damage_band_begin, damage_band_end);
      do_supg_computation = !hot_points.empty();
    }
    if (do_supg_computation)
    {
      xIntegrationRule& damage_band_integration_rule = *geom.getIntegrationRule(Geom::DAMAGEBAND, integration_order);

      const xMesh* front_mesh = geom.getFrontMeshIn();
Kevin Moreau's avatar
Kevin Moreau committed
269 270
      xExportGmshAscii pexport;
      Export(front_mesh, pexport, "front_mesh");
271 272
#ifdef MODEBND
      const xMesh* geo_bnd_mesh = geom.getFrontMeshGeo();
Kevin Moreau's avatar
Kevin Moreau committed
273 274
      Export(geo_bnd_mesh, pexport, "geo_bnd_mesh");
      Export(ls_in.getSupport().getMesh(), pexport, "mesh");
275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
#endif

#ifdef VECTORLS
      xUpperCreatorFiltered upper_creator_filtered(geom.getFilter(Geom::DAMAGEBANDSUPPORT));
      xCreator upper_creator_front_filtered;
#else
      xUpperCreatorFilteredRecursive upper_creator_filtered(geom.getFilter(Geom::DAMAGEBANDSUPPORT));
      xUpperCreatorFilteredRecursive upper_creator_front_filtered(geom.getFilter(Geom::FRONTINSUPPORT));
      // xUpperCreatorFilteredRecursive upper_creator_front_filtered(geom.getFilter(Geom::FRONTINSUPPORTELEMENTWISE));
#endif
#ifdef MODEBND
      // xEntityFilter damage_band_support_filter = geom.getFilter(Geom::DAMAGEBANDSUPPORT);
      xEvalLevelSet<xIdentity<double> > eval_ls_in(ls_in);
      xEvalGradSmoothLevelSet<xIdentity<xVector> > eval_grad_ls_geo(*ls_geo);
      xEvalGradSmoothLevelSet<xIdentity<xVector> > eval_grad_ls_in(ls_in);
      std::vector<mEntity*> front_mesh_region;
      const double zero = 0.;
      for(xIter it=front_mesh->begin(1); it!=front_mesh->end(1); ++it) {
        front_mesh_region.push_back(*it);
      }

      std::vector<mEntity*> front_and_bnd_mesh_region(front_mesh_region);
297
      xSubMesh& bnd_sub_mesh = geo_bnd_mesh->createSubMesh("bnd");
298 299 300 301 302 303 304 305 306 307 308
      xCreator creator;
      for(xIter it=geo_bnd_mesh->begin(1); it!=geo_bnd_mesh->end(1); ++it) {
        mEntity* e = *it;
        mEntity* ec = creator(e);
        xGeomElem geom_elem(ec);
        bool is_damaged = true;
        for(int i=0; i<e->size(0) && is_damaged; ++i) {
          mVertex* v = static_cast<mVertex*>(e->get(0,i));
          geom_elem.setUVWForXYZ(v->point());
          double val;
          eval_ls_in(&geom_elem, &geom_elem, val);
309
          is_damaged = is_damaged && (val>=0.) && (val<=lc);
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
        }        
        if (is_damaged) {
          xPartition partition;
          geo_bnd_mesh->getPartition(e, partition);
          for(xPartition::iterator itp=partition.begin(); itp!=partition.end(); ++itp) {
            mEntity* ep = *itp;
            double x=0;
            double y=0;
            for(int i=0; i<ep->size(0); ++i) {
              mVertex* v = static_cast<mVertex*>(ep->get(0,i));
              x+=v->point()(0);
              y+=v->point()(1);
            }
            mPoint p(x/static_cast<double>(ep->size(0)), y/static_cast<double>(ep->size(0)), 0.);
            geom_elem.setUVWForXYZ(p);
            xVector vect, normal;
            eval_grad_ls_in(&geom_elem, &geom_elem, vect);
            eval_grad_ls_geo(&geom_elem, &geom_elem, normal);
            if(vect*normal<zero) {
329
              // front_and_bnd_mesh_region.push_back(ep);
330 331 332 333 334
              bnd_sub_mesh.add(ep);
            }
          }
        }
      }
335 336 337 338 339
      if(bnd_sub_mesh.size(1)>1) {
        for (xIter it=bnd_sub_mesh.begin(1); it!=bnd_sub_mesh.end(1); ++it) {
          front_and_bnd_mesh_region.push_back(*it);
        }
      }
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
      if(bnd_sub_mesh.size(1)) {
        std::ostringstream oss;
        oss.width(5);
        oss.fill('0');
        oss << step_;
        bnd_sub_mesh.modifyAllState();
        bnd_sub_mesh.exportGmsh("bnd_"+oss.str()+".msh");
      }
      geo_bnd_mesh->deleteSubMesh("bnd");
#else
#ifndef VECTORLS
      AcceptEntityToEntity filter(upper_creator_filtered, geom.getFilter(Geom::DAMAGEBANDSUPPORT));
      xFilteredRegion<xIter, xEntityFilter> front_mesh_region(front_mesh->begin(dim-1), front_mesh->end(dim-1), filter);
#endif
#endif

      SolverNonLocal solver_nonlocal;
      lalg::xDenseMatrix P;
      xEvalConstant<double> eval_supg_stabilization(supg_stabilization_factor*min_element_size);
      xIntegrationRuleBasic integration_rule_basic(2);
      front_mode_extension.computeModesFromFrontModesWithSUPG(eval_normalized_grad_ls, eval_supg_stabilization,
                                                              solver_nonlocal,
                                                              damage_band_integration_rule, integration_rule_basic,
                                                              damage_band_begin, damage_band_end,
#ifdef MODEBND
                                                              front_and_bnd_mesh_region.begin(), front_and_bnd_mesh_region.end(),
#else
#ifdef VECTORLS
                                                              geom.beginFront(), geom.endFront(),
#else
                                                              front_mesh_region.begin(), front_mesh_region.end(),
#endif
#endif
                                                              upper_creator_filtered, upper_creator_front_filtered,
                                                              mode_min_support_length_factor*min_element_size,
                                                              P, step_);

      lalg::xCSRVector front_velocity(0);
#ifndef ESHELBY
      double nonlocal_dissipation = eval_front_velocity.computeLumped(eval_energy_release_rate_, eval_damage_derivative,
                                                                      eval_damage_rate_wrapper_,
                                                                      P,
                                                                      damage_band_integration_rule,
                                                                      damage_band_begin, damage_band_end,
                                                                      front_velocity);
#else
      lalg::xDenseMatrix Q;
      front_mode_extension.computeVectorModesFromModes(integration_rule_basic,
                                                       front_mesh_region.begin(), front_mesh_region.end(),
                                                       upper_creator_front_filtered,
                                                       P, Q);

      const xMesh* crack_mesh = geom.getFrontMeshOut();
      double nonlocal_dissipation = eval_front_velocity.computeLumped(eval_eshelby_tensor, eval_damage_derivative, eval_normalized_grad_ls,
                                                                      eval_damage_rate_wrapper_,
                                                                      P, Q,
                                                                      damage_band_integration_rule,
                                                                      damage_band_begin, damage_band_end,
                                                                      front_mesh_region.begin(), front_mesh_region.end(),
                                                                      crack_mesh->begin(1), crack_mesh->end(1),
                                                                      upper_creator_front_filtered,
                                                                      front_velocity);
#endif
      cumulate_nonlocal_dissipation+=time_step_*nonlocal_dissipation;
      double max_val = *std::max_element(front_velocity.begin(), front_velocity.end());
      to_move = max_val>0.;
      post_processing.exportOnSpace("front_velocity", step_, eval_front_velocity, damage_band_integration_rule, damage_band_begin, damage_band_end);
      post_processing.exportOnTime("max_front_velocity", time_, max_val);
      post_processing.exportOnTime("nonlocal_dissipation", time_, nonlocal_dissipation);
      post_processing.exportOnTime("cumulate_nonlocal_dissipation", time_, cumulate_nonlocal_dissipation);
      observer.setToDo("reassemble_mass", to_move);

      if (to_move)
      {
414 415 416 417 418 419 420
#ifdef KALTHOFF_HARD_CODED_TEST
        std::map<int, double> max_val_map;
        sensors.measureMaxOnAllBoxes(eval_front_velocity, damage_band_integration_rule, damage_band_begin, damage_band_end, max_val_map);
        post_processing.exportOnTime("max_front_velocity_box1", time_, max_val_map[1]);
        post_processing.exportOnTime("max_front_velocity_box2", time_, max_val_map[2]);
#endif

421 422
        narrow_band = &(ls_in.getSupport().getMesh()->createSubMesh("narrow_band"));
        BuildBand(geom.begin(Geom::LINEAR, dim), geom.end(Geom::LINEAR, dim), -lc, 0., ls_in, *narrow_band, narrow_band_attach_tag);
423

424 425
        // xAcceptOnIntTag filter_narrow_band_support(narrow_tag, 1);
        xAcceptInSubMesh filter_narrow_band_support(*narrow_band);
426 427 428 429 430 431 432 433
        xUpperCreatorFilteredRecursive upper_creator_narrow_band_filtered(filter_narrow_band_support);
        xEntityFilter filter_narrow_band = geom.getFilter(Geom::FRONTCOVERSIN); 
        xIntegrationRulePartition narrow_band_integration_rule(filter_narrow_band, integration_order);
#ifdef MODEBND
        eval_front_velocity_extended.compute(eval_minus_normalized_grad_ls, eval_supg_stabilization,
                                             eval_front_velocity,
                                             solver_nonlocal,
                                             narrow_band_integration_rule, integration_rule_basic,
434
                                             narrow_band->begin(), narrow_band->end(),
435 436 437 438 439 440 441 442
                                             front_mesh_region.begin(), front_mesh_region.end(),
                                             upper_creator_narrow_band_filtered, upper_creator_front_filtered,
                                             mode_min_support_length_factor*min_element_size);
#else
        lalg::xDenseMatrix P2;
        front_mode_extension.computeModesFromFrontModesWithSUPG(eval_minus_normalized_grad_ls, eval_supg_stabilization,
                                                                solver_nonlocal,
                                                                narrow_band_integration_rule, integration_rule_basic,
443
                                                                narrow_band->begin(), narrow_band->end(),
444 445 446 447 448 449 450 451 452 453 454
#ifdef VECTORLS
                                                                geom.beginFront(), geom.endFront(),
#else
                                                                front_mesh_region.begin(), front_mesh_region.end(),
#endif
                                                                upper_creator_narrow_band_filtered, upper_creator_front_filtered,
                                                                mode_min_support_length_factor*min_element_size,
                                                                P2, step_);

        eval_front_velocity_extended.compute(P2,
                                             front_velocity,
455
                                             narrow_band->begin(), narrow_band->end());
456 457
#endif

458
        post_processing.exportOnSpace("front_velocity_extended", step_, eval_front_velocity_extended, narrow_band_integration_rule, narrow_band->begin(), narrow_band->end());
459
      }
460 461 462 463 464
      else {
#ifdef KALTHOFF_HARD_CODED_TEST
        post_processing.exportOnTime("max_front_velocity_box1", time_, 0.);
        post_processing.exportOnTime("max_front_velocity_box2", time_, 0.);
#endif
465 466 467 468 469 470 471 472 473 474 475 476
      }
    }
  }
}

void TLSSolverMode::move(double time_, int step_, double time_step_)
{
  if (to_move)
  {
    xScale<double> time_step_scale(time_step_);
    xEvalUnary<xScale<double> > eval_delta_ls_in_extended(time_step_scale, eval_front_velocity_extended);
    xEvalUnary<xScale<double> > eval_delta_ls_in(time_step_scale, eval_front_velocity);
477 478
    xSubMesh* damage_band = &(ls_in.getSupport().getMesh()->createSubMesh("damage_band"));
    BuildBand(geom.begin(Geom::DAMAGEBANDSUPPORT, dim), geom.end(Geom::DAMAGEBANDSUPPORT, dim), 0., lc, ls_in, *damage_band, damage_band_attach_tag);
479 480 481 482 483 484 485 486
    PartialShiftLevelSetModifier narrow_band_shifter(eval_delta_ls_in_extended, *narrow_band, narrow_band_attach_tag);
    PartialShiftLevelSetModifier damage_band_shifter(eval_delta_ls_in, *damage_band, damage_band_attach_tag);
    ls_in.accept(narrow_band_shifter);
    ls_in.accept(damage_band_shifter);
    ClearAttachedEntity(narrow_band->begin(0), narrow_band->end(0), narrow_band_attach_tag);
    ClearAttachedEntity(damage_band->begin(0), damage_band->end(0), damage_band_attach_tag);
    RemoveIsolatedNegativeLevelSet<xIter> ls_remove_isolated(ls_in.getSupport().begin(2), ls_in.getSupport().end(2));
    ls_in.accept(ls_remove_isolated);
487 488 489 490 491 492

    // // tmp
    // std::ostringstream oss;
    // oss.width(5);
    // oss.fill('0');
    // oss << step_;
493 494 495 496 497 498 499 500 501 502 503
    // std::ofstream outf("damage_band_"+oss.str()+".geo", std::ofstream::out);
    // for (xIter it=damage_band->begin(0); it!=damage_band->end(0); ++it) {
    //   mPoint p = static_cast<mVertex*>(*it)->point();
    //   outf << "Point ( "+toString((*it)->getId())+" ) = { "+toString(p(0))+" , "+toString(p(1))+" , "+toString(p(2))+" };" << std::endl;
    // }
    // std::ofstream outf2("narrow_band_"+oss.str()+".geo", std::ofstream::out);
    // for (xIter it=narrow_band->begin(0); it!=narrow_band->end(0); ++it) {
    //   mPoint p = static_cast<mVertex*>(*it)->point();
    //   outf2 << "Point ( "+toString((*it)->getId())+" ) = { "+toString(p(0))+" , "+toString(p(1))+" , "+toString(p(2))+" };" << std::endl;
    // }
    // damage_band->modifyAllState();
504
    // damage_band->exportGmsh("damage_band_"+oss.str()+".msh");
505 506
    // narrow_band->modifyAllState();
    // narrow_band->exportGmsh("narrow_band_"+oss.str()+".msh");
507 508 509 510 511 512
    // // tmp

    ls_in.getSupport().getMesh()->deleteSubMesh("narrow_band");
    ls_in.getSupport().getMesh()->deleteSubMesh("damage_band");
    narrow_band = NULL;
    damage_band = NULL;
513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540
#ifdef MODEBND
    updateDomain();
    xEntityFilter damage_band_support_filter = geom.getFilter(Geom::DAMAGEBANDSUPPORT);
    const xMesh* front_mesh = geom.getFrontMeshIn();
    const xMesh* geo_bnd_mesh = geom.getFrontMeshGeo();
    xEvalLevelSet<xIdentity<double> > eval_ls_in(ls_in);
    xEvalGradSmoothLevelSet<xIdentity<xVector> > eval_grad_ls_geo(*ls_geo);
    xEvalGradSmoothLevelSet<xIdentity<xVector> > eval_grad_ls_in(ls_in);
    std::set<mEntity*> container;
    const double zero = 0.;
    for(xIter it=front_mesh->begin(0); it!=front_mesh->end(0); ++it) {
      container.insert(*it);
    }
    xSubMesh& bnd_sub_mesh = geo_bnd_mesh->createSubMesh("bnd");
    xCreator creator;
    for(xIter it=geo_bnd_mesh->begin(1); it!=geo_bnd_mesh->end(1); ++it) {
      mEntity* e = *it;
      mEntity* ec = creator(e);
      xGeomElem geom_elem(ec);
      xPartition partition;
      geo_bnd_mesh->getPartition(e, partition);
      for(xPartition::iterator itp=partition.begin(); itp!=partition.end(); ++itp) {
        mEntity* ep = *itp;
        for(int i=0; i<ep->size(0); ++i) {
          mVertex* v = static_cast<mVertex*>(ep->get(0,i));
          geom_elem.setUVWForXYZ(v->point());
          double val;
          eval_ls_in(&geom_elem, &geom_elem, val);
541 542 543 544 545 546 547 548
          if(val>=0. && val<=lc) {
            // xVector vect, normal;
            // eval_grad_ls_in(&geom_elem, &geom_elem, vect);
            // eval_grad_ls_geo(&geom_elem, &geom_elem, normal);
            // if(vect*normal<zero) {
            container.insert(v);
            bnd_sub_mesh.add(v);
            // }
549 550 551 552 553 554 555 556 557 558 559 560 561 562
          }
        }
      } 
    }
    if(bnd_sub_mesh.size(0)) {
      std::ostringstream oss;
      oss.width(5);
      oss.fill('0');
      oss << step_;
      bnd_sub_mesh.modifyAllState();
      bnd_sub_mesh.exportGmsh("bnd_reinit_"+oss.str()+".msh");
    }
    geo_bnd_mesh->deleteSubMesh("bnd");
    xtls::reInitLevelSetFromFrontAndBoundaryBruteForce ls_reinit2;
563
    SetTagVertices(geom.begin(Geom::DAMAGED,2), geom.end(Geom::DAMAGED,2), 111112);
564
    xAcceptOnIntTag damage_band_filter(111112, 1);
565
    xUpperCreatorFilteredRecursive upper_creator_filtered_recursive(damage_band_support_filter);
566
    ls_reinit2.setBoundaryValues(eval_ls_in, container.begin(), container.end(),
567
                                 upper_creator_filtered_recursive, damage_band_filter, xAcceptAll()); // geom.getFilter(Geom::FRONTINSUPPORT));
568
    ls_reinit2.visit(ls_in, ls_in.getSupport());
569
    ClearTagVertices(geom.begin(Geom::DAMAGED,2), geom.end(Geom::DAMAGED,2), 111112);
570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
    xExportGmshAscii pexport;
    Export(ls_in, pexport, "ls_in_reinit");
    container.clear();
#else
    ls_in.accept(ls_reinit);
#endif
    updateDomain();
    geom.incrementFrontAdvance();
    to_move=false;
  }
}

void TLSSolverMode::updateDomain()
{
  geom.clear();
  try {
    geom.updateDomain(ls_in.getSupport().getMesh(), ls_in, lc, ls_geo);
  }
588 589
  catch(...) {
    std::cout << "Warning: catch error numbered " << std::endl;
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640
    smoothLevelSet();
    geom.clear();
    geom.updateDomain(ls_in.getSupport().getMesh(), ls_in, lc, ls_geo);
  }
#ifdef VECTORLS
  eval_adim_ls.update(boost::bind(geom.getPartionFunction(), _1, _2, geom.getFilter(Geom::DAMAGEBAND)), geom.getFrontMeshOut());
#endif
  eval_damage.setFilters(geom.getFilter(Geom::DAMAGEBAND), geom.getFilter(Geom::SANEMATTER));
  eval_damage_derivative.setFilters(geom.getFilter(Geom::DAMAGEBAND), geom.getFilter(Geom::SANEMATTER));
  eval_damage_rate.setFilter(geom.getFilter(Geom::DAMAGEBAND));
}

void TLSSolverMode::smoothLevelSet()
{
  xLevelSetSmoother ls_smoother;
  xLevelSet ls_copy(ls_in);
  ls_smoother.visit(ls_copy,ls_in);
}

void TLSSolverMode::exportGmsh(double time_, int step_, double time_step_, const IterDomain& begin_, const IterDomain& end_, const xIntegrationRule& integration_rule_)
{
  if (search_for_initiation.isFirstInitiated())
  {
    post_processing.exportOnTime("damage", time_, eval_damage);
    post_processing.exportOnSpace("damage", step_, eval_damage, integration_rule_, begin_, end_);
    post_processing.exportOnSpace("ls_in", step_, ls_in);
    xEvalGradLevelSet<xIdentity<xVector> > eval_grad_ls_in(ls_in);
    post_processing.exportOnSpace("ls_in_grad", step_, eval_grad_ls_in, integration_rule_, begin_, end_);
  }
  post_processing.storeLevelSet("ls_in", step_, ls_in);
  // post_processing.storeValue("cumulate_dissipation",          step_, cumulate_dissipation);
  // post_processing.storeValue("cumulate_fake_dissipation",     step_, cumulate_fake_dissipation);
  post_processing.storeValue("cumulate_nonlocal_dissipation", step_, cumulate_nonlocal_dissipation);
  post_processing.storeValue("newly_initiated",               step_, search_for_initiation.isNewlyInitiated());
  post_processing.storeValue("first_initiated",               step_, search_for_initiation.isFirstInitiated());
}

const xEval<double>& TLSSolverMode::getEvalDamage()
{ 
  return eval_damage;
}

const xEval<double>& TLSSolverMode::getEvalDamageRate()
{
  return eval_damage_rate;
}

const xEval<xVector>& TLSSolverMode::getEvalGradLs()
{ 
  return eval_normalized_grad_ls;
}