/test/api/egg_ut.cpp

//
// avro - Adaptive Voronoi Remesher
//
// Copyright 2017-2022, Philip Claude Caplan
// All rights reserved
//
// Licensed under The GNU Lesser General Public License, version 2.1
// See http://www.opensource.org/licenses/lgpl-2.1.php
//
#include "unit_tester.hpp"

#include "geometry/egads/context.h"
#include "geometry/egads/data.h"
#include "geometry/egads/object.h"

#include "graphics/application.h"
#include "graphics/math.h"

#include "library/ckf.h"
#include "numerics/geometry.h"
#include "numerics/linear_algebra.h"

#include "avro.h"
#include "avro_config.h"

#include <egads.h>

#include <cmath>

using namespace avro;

UT_TEST_SUITE(api_custom_geometry_test_suite)

#if AVRO_NO_ESP == 0

void
add_children( EGADSGeneralGeometry& egg , Entity* entity , std::map<ego,ego>& duplicates ) {

  ego e = static_cast<EGADS::Object*>(entity)->object();

  // do not add the duplicate ego's into the geometry hierarchy
  if (duplicates.find(e) != duplicates.end()) return;

  // add the children of this entity
  for (index_t k = 0; k < entity->nb_children(); k++) {
    ego ek = static_cast<EGADS::Object*>(&entity->child(k))->object();

    // again, do not add the duplicate ego's into the geometry hierarchy
    if (duplicates.find(ek) != duplicates.end()) continue;

    // add the child and recursively add its children as well
    egg.add_child(e,ek);
    add_children(egg,&entity->child(k),duplicates);
  }
}

UT_TEST_CASE(test1)
{
  // this example glues two boxes together and tags the interior Face between them as an "interior"
  // (note: interior Edges and Nodes on this Face are also tagged as interior)
  // please follow the different steps outlined below to set up your own custom geometry
  // (see the STEP ... and END STEP)

  // you can pass in an existing EGADS context if you want
  EGADS::Context ctx;

  // STEP 1: initialize  the geometry and avro context
  const coord_t dim = 3;
  coord_t number = dim;
  coord_t udim = dim -1;
  EGADSGeneralGeometry egg(ctx.get(),number-1);
  avro::Context adapter(dim,dim,udim,false);
  // END STEP 1

  // create the first box
  ego box1;
  real_t params1[6] = { 0.0 , 0.0 , 0.0 , 1.0 , 1.0 , 0.5 };
  EGADS_ENSURE_SUCCESS( EG_makeSolidBody( egg.context()->get() , BOX , params1 , &box1 ) );

  // create the second box
  ego box2;
  real_t params2[6] = { 0.0 , 0.0 , 0.5 , 1.0 , 1.0 , 0.5 };
  EGADS_ENSURE_SUCCESS( EG_makeSolidBody( egg.context()->get() , BOX , params2 , &box2 ) );

  // create the avro bodies (this is only needed for this demo, ideally you would have your own wrapper around ego's)
  avro::EGADS::Body b1( *egg.context() , box1 );
  avro::EGADS::Body b2( *egg.context() , box2 );
  b1.build_hierarchy();
  b2.build_hierarchy();

  // add each body ego to the egg
  egg.add_body( box1 );
  egg.add_body( box2 );

  // retrieve the entities of the first body
  std::vector<Entity*> entities1;
  b1.get_entities(entities1);
  printf("body1 has %lu entities\n",entities1.size());

  // retrieve the entities of the second body
  std::vector<Entity*> entities2;
  b2.get_entities(entities2);
  printf("body2 has %lu entities\n",entities2.size());

  // STEP 2: add the ego's to the geometry
  // add the ego's for each body (even if some are duplicates and won't be used)
  std::vector<Entity*> entities;
  for (index_t k = 0; k < entities1.size(); k++) {
    ego ek = static_cast<EGADS::Object*>(entities1[k])->object();
    egg.add_object( box1 , ek );
    entities.push_back( entities1[k] );
  }

  for (index_t k = 0; k < entities2.size(); k++) {
    ego ek = static_cast<EGADS::Object*>(entities2[k])->object();
    egg.add_object( box2 , ek );
    entities.push_back( entities2[k] );
  }
  printf("--> added objects\n");
  // END STEP 2

  // generate a mesh in each box
  index_t N = 5;
  CKF_Triangulation mesh1( {N,N,N} );
  CKF_Triangulation mesh2( {N,N,N} );

  // squish the z-coordinates
  avro_assert( mesh1.points().nb() == mesh2.points().nb() );
  for (index_t k = 0; k < mesh1.points().nb(); k++) {
    mesh1.points()[k][2] *= 0.5;
    mesh2.points()[k][2] = mesh1.points()[k][2] + 0.5;
  }

  // attach the mesh points to the bodies (again, this is only for the demo)
  // if you created your own meshes to start
  // (say from TetGen, then you should retain which ego each vertex is on)
  mesh1.points().attach( b1 );
  mesh2.points().attach( b2 );

  // map duplicate entities (which should be at z = 0.5)
  // there should be 1 face, 4 edges and 4 nodes which are duplicates
  // (see the topological number printed below)
  std::map<Entity*,Entity*> entity_map;
  std::map<ego,ego> ego_map;
  for (index_t i = 0; i < entities1.size(); i++)
  for (index_t j = 0; j < entities2.size(); j++) {

    ego ei = static_cast<EGADS::Object*>(entities1[i])->object();
    ego ej = static_cast<EGADS::Object*>(entities2[j])->object();

    int same = EG_isSame(ei,ej);
    if (same == EGADS_SUCCESS)
      printf("number = %u, i = %lu, j = %lu are the same!\n",entities1[i]->number(),i,j);
    else
      continue;

    // we will only keep entities1[i]
    entity_map.insert( {entities2[j],entities1[i]} );
    ego_map.insert( {ej,ei} );

    // mark the entity as interior
    entities1[i]->set_interior(true);

    // add the unique entity as a child of the duplicate entity's parents
    std::shared_ptr<Entity> ep = b1.lookup(ei); // retrieve the entity that we want to keep
    avro_assert( ep != nullptr );
    for (index_t k = 0; k < entities2[j]->nb_parents(); k++) {
      entities2[j]->parents(k)->add_child(ep);
    }
  }
  printf("--> number of duplicate entities = %lu\n",entity_map.size());

  // STEP 3: set the interior entities for the egg
  for (index_t k = 0; k < entities1.size(); k++) {
    if (!entities1[k]->interior()) continue;
    ego e = static_cast<EGADS::Object*>(entities1[k])->object();
    egg.set_interior(e);
  }
  // END STEP 3

  // STEP 4: add the ego parent-child relationships
  // add the adjacency relationships
  for (index_t k = 0; k < entities.size(); k++) {
    add_children(egg,entities[k],ego_map);
  }
  printf("--> adjacencies added\n");
  // END STEP 4

  // STEP 5: finalize the customized geometry and pass it to the context
  egg.finalize();
  adapter.define_geometry(egg);
  // END STEP 5

  // STEP 6: generate your mesh from the geometry
  // (in general you might use TetGen and retain the ego's each vertex is on)
  // create one final mesh
  Points points(dim);
  Topology<Simplex> topology(points,number);

  // copy the points
  mesh1.points().copy( points );
  std::map<index_t,index_t> point_map;
  for (index_t k = 0; k < mesh2.points().nb(); k++) {

    Entity* ek = mesh2.points().entity(k);

    // check if the point is not a duplicate
    if (ek == nullptr || entity_map.find(ek) == entity_map.end()) {
      index_t idx = points.nb();
      point_map.insert( {k,idx} );
      points.create( mesh2.points()[k] );
      points.set_entity( idx , ek );
      points.set_param( idx , mesh2.points().u(k) );
      continue;
    }

    // determine which point in mesh1 this is closest to
    int idx = -1;
    for (index_t i = 0; i < mesh1.points().nb(); i++) {

      if (mesh1.points().entity(i) == nullptr) continue;
      if (!mesh1.points().entity(i)->interior()) continue; // reduce the loop down to O(# duplicate points)

      real_t d = numerics::distance( mesh1.points()[i] , mesh2.points()[k] , dim );
      if (d < 1e-12) {
        idx = i;
        break;
      }
    }
    avro_assert( idx >= 0 );
    point_map.insert( {k,idx} );
  }
  printf("--> number of duplicate points = %lu\n",
                     mesh1.points().nb() + mesh2.points().nb() - points.nb());
  UT_ASSERT_EQUALS( mesh1.points().nb() + mesh2.points().nb() - points.nb() , N*N );

  // compute the parameter coordinates (which may have changed when we adjusted for duplicates)
  for (index_t k = 0; k < points.nb(); k++) {
    Entity* e = points.entity(k);
    if (e == nullptr) continue;
    std::vector<real_t> X( points[k] , points[k] + dim );
    std::vector<real_t> U( dim-1 );
    e->inverse( X , U );
    points.set_param( k , U.data() );
  }
  // END STEP 6

  // check the number of points on geometry Nodes is 12 (8 + 8 - the 4 duplicates)
  index_t nb_nodes = 0;
  for (index_t k = 0; k < points.nb(); k++) {
    if (points.entity(k) == nullptr) continue;
    if (points.entity(k)->number() == 0) nb_nodes++;
  }
  UT_ASSERT_EQUALS( nb_nodes , 12 );

  // copy the simplices
  for (index_t k = 0; k < mesh1.nb(); k++)
    topology.add( mesh1(k) , mesh1.nv(k) );
  for (index_t k = 0; k < mesh2.nb(); k++) {
    std::vector<index_t> elem( mesh2(k) , mesh2(k) + mesh2.nv(k) );
    for (index_t i = 0; i < elem.size(); i++)
      elem[i] = point_map.at( elem[i] );
    topology.add( elem.data() , elem.size() );
  }

  // STEP 7: retrieve the association between ego's and the context's integer labels
  std::map<ego,int> ids;
  adapter.get_ego_ids(ids);
  printf("--> nb geometry ids = %lu\n",ids.size());

  // determine the geometry (integer) label for every vertex
  std::vector<real_t> coordinates( dim*points.nb() );
  std::vector<int> geometry( points.nb() , -1 );
  std::vector<real_t> parameters( points.nb() * (dim-1) , unset_value );
  for (index_t k = 0; k < points.nb(); k++) {
    for (coord_t d = 0; d < dim; d++)
      coordinates[k*dim+d] = points[k][d];
    Entity* entity = points.entity(k);
    if (entity == nullptr) continue;
    parameters[2*k]   = points.u(k)[0];
    parameters[2*k+1] = points.u(k)[1];
    ego e = static_cast<EGADS::Object*>(entity)->object();
    geometry[k] = ids.at(e);
  }
  // END STEP 7, then proceed to use the context as usual

  // send the mesh and geometry information to the context
  adapter.load_mesh( coordinates , topology.data() );
  adapter.load_geometry( geometry , parameters );

  index_t nb_points = points.nb();
  index_t nb_rank = number*(number+1)/2;
  std::vector<real_t> metrics( nb_points * nb_rank , 0.0 );

  // constant diagonal metric requesting a size of h
  real_t h = 0.125;
  for (index_t k = 0; k < nb_points; k++) {
    index_t idx = k*nb_rank;
    for (index_t i = 0; i < number; i++) {
      metrics[idx] = 1./(h*h);
      idx += number - i;
    }
  }

  //adapter.parameters().set( "output redirect" , "avro-output.txt" );
  adapter.parameters().set( "curved" , false );

  // perform the adaptation
  adapter.adapt(metrics);

  // retrieve the mesh from the context
  std::vector<index_t> elements;
  adapter.retrieve_mesh( coordinates , elements );
  nb_points = coordinates.size() / dim;
  index_t nb_elements = elements.size() / (number+1);
  printf("adapted mesh has %lu points and %lu elements\n",nb_points,nb_elements);

  // retrieve the geometry metadata
  std::vector<real_t> param;
  std::vector<int> entity_idx;
  adapter.retrieve_geometry( entity_idx , param );
  assert( param.size() == udim*nb_points );
  assert( entity_idx.size() == nb_points );
  for (index_t k = 0; k < nb_points; k++) {
    printf("\tpoint [%lu] (%g,%g,%g) on entity %d: u = ( ",k,coordinates[3*k],coordinates[3*k+1],coordinates[3*k+2],entity_idx[k]);
    for (index_t j = 0; j < udim; j++)
      printf("%g ",param[k*udim+j]);
    printf(")\n");
  }

  // retrieve the boundary facets
  std::vector<std::vector<index_t>> facets;

  // either of the following two functions will work (the first uses the current mesh stored in the context, the second uses the geometry information of each vertex along with the mesh elements to infer the boundary)
  // note that vertex coordinates are not required since this is a purely topological operation, which is why the second method does not require vertex coordinates
  //context.retrieve_boundary( facets , geometry );
  adapter.retrieve_boundary( entity_idx , elements , facets , geometry , true ); // true because we have interior entities

  assert( facets.size() == geometry.size() );
  std::vector<const real_t*> x(number);
  real_t area = 0.0;
  for (index_t bc = 0; bc < facets.size(); bc++) {
    const std::vector<index_t>& facets_on_bc = facets[bc];
    index_t nb_facets = facets_on_bc.size() / number;
    printf("there are %lu facets on bc %lu with geometry id %d\n",nb_facets,bc,geometry[bc]);

    for (index_t k = 0; k < nb_facets; k++) {


      // calculate the area of the facet
      for (index_t j = 0; j < number; j++)
        x[j] = &coordinates[3*facets_on_bc[k*number+j]];

      area += numerics::volume_nd( x , dim );
    }
  }

  // we have 2 external 1x1 faces, 8 external 0.5x1 faces, 1 internal 1x1 face
  printf("area = %g\n",area);
  UT_ASSERT_NEAR( area , 7.0 , 1e-12 );

}
UT_TEST_CASE_END(test1)

#endif // AVRO_NO_ESP

UT_TEST_SUITE_END(api_custom_geometry_test_suite)