Skip to content

Commit

Permalink
lsMesh/lsPointData as templates (#37)
Browse files Browse the repository at this point in the history 
* Changed lsPointData and lsMesh to templates, so they can be used for float or double.

* Fixed numerical issue in lsMakeGeometry<float> when creating spheres. 

* Added warning output for null vectors in lsCalculateNormalVectors.

* Removed grid spacing array from lsToDiskMesh.
  • Loading branch information
@XaverKlemenschits
XaverKlemenschits authored Mar 8, 2021
1 parent 1e9953f commit f58b429
Show file tree
Hide file tree
Showing 276 changed files with 3,108 additions and 2,669 deletions.
11 changes: 8 additions & 3 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,11 @@ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC")
endif()

# set whether to build static versions
option(VIENNALS_STATIC_BUILD "Build all targets with static links." OFF)
option(VIENNALS_STATIC_BUILD "Link dependencies statically." OFF)
if(VIENNALS_STATIC_BUILD)
set(VTK_DIR $ENV{VTK_STATIC_DIR})
if(NOT VTK_DIR)
set(VTK_DIR $ENV{VTK_STATIC_DIR})
endif(NOT VTK_DIR)
list(APPEND VIENNALS_LIBRARIES "-static")
endif(VIENNALS_STATIC_BUILD)

Expand All @@ -50,7 +52,10 @@ if(VIENNALS_USE_VTK)
include(${VTK_USE_FILE})
# only link needed vtk libraries for static builds
if(VIENNALS_STATIC_BUILD)
set(VTK_LIBRARIES vtksys;vtkIOCore;vtkexpat;vtklz4;vtkzlib;vtklzma;vtkdoubleconversion;vtkCommonMisc;vtkCommonSystem;vtkIOXML)
set(VTK_LIBRARIES vtksys;vtkIOCore;vtkexpat;vtklz4;vtkzlib;vtklzma;vtkdoubleconversion;vtkCommonMisc;vtkCommonSystem;vtkIOXML;vtkFiltersCore;vtkFiltersGeneral;vtkFiltersGeometry;vtkCommonExecutionModel;vtkCommonDataModel;vtkCommonColor;vtkIOXMLParser;vtkIOCore;vtkCommonMisc;vtkCommonComputationalGeometry;vtkCommonTransforms;vtkCommonMath;)

message(STATUS "Static build: Only linking ${VTK_LIBRARIES}")
# vtksys;vtkCommonCore;vtkCommonMath;vtkCommonMisc;vtkCommonSystem;vtkCommonTransforms;vtkCommonDataModel;vtkCommonColor;vtkCommonExecutionModel;vtkCommonComputationalGeometry;vtkFiltersCore;vtkFiltersGeneral;vtkImagingCore;vtkImagingFourier;vtkFiltersStatistics;vtkFiltersExtraction;vtkInfovisCore;vtkFiltersGeometry;vtkFiltersSources;vtkRenderingCore;vtkzlib;vtkfreetype;vtkRenderingFreeType;vtkRenderingContext2D;vtkChartsCore;vtkDICOMParser;vtkdoubleconversion;vtklz4;vtklzma;vtkIOCore;vtkIOLegacy;vtkexpat;vtkIOXMLParser;vtkDomainsChemistry;vtkglew;vtkRenderingOpenGL2;vtkDomainsChemistryOpenGL2;vtkIOXML;vtkParallelCore;vtkFiltersAMR;vtkFiltersFlowPaths;vtkFiltersGeneric;vtkImagingSources;vtkFiltersHybrid;vtkFiltersHyperTree;vtkImagingGeneral;vtkFiltersImaging;vtkFiltersModeling;vtkFiltersParallel;vtkFiltersParallelImaging;vtkFiltersPoints;vtkFiltersProgrammable;vtkFiltersSMP;vtkFiltersSelection;vtkFiltersTexture;vtkFiltersTopology;verdict;vtkFiltersVerdict;vtkmetaio;vtkjpeg;vtkpng;vtktiff;vtkIOImage;vtkImagingHybrid;vtkInfovisLayout;vtkInteractionStyle;vtkImagingColor;vtkRenderingAnnotation;vtkRenderingVolume;vtkInteractionWidgets;vtkViewsCore;vtklibproj;vtkGeovisCore;vtkhdf5_src;vtkhdf5_hl_src;vtkIOAMR;vtkIOAsynchronous;vtkpugixml;vtkIOCityGML;vtkIOEnSight;vtknetcdf;vtkexodusII;vtkIOExodus;vtkgl2ps;vtkRenderingGL2PSOpenGL2;vtkIOExport;vtkIOExportOpenGL2;vtklibharu;vtkIOExportPDF;vtkIOGeometry;vtkIOImport;vtklibxml2;vtkIOInfovis;vtkIOLSDyna;vtkIOMINC;vtkogg;vtktheora;vtkIOMovie;vtkIONetCDF;vtkIOPLY;vtkjsoncpp;vtkIOParallel;vtkIOParallelXML;vtksqlite;vtkIOSQL;vtkIOSegY;vtkIOTecplotTable;vtkIOVeraOut;vtkIOVideo;vtkImagingMath;vtkImagingMorphological;vtkImagingStatistics;vtkImagingStencil;vtkInteractionImage;vtkRenderingContextOpenGL2;vtkRenderingImage;vtkRenderingLOD;vtkRenderingLabel;vtkRenderingVolumeOpenGL2;vtkViewsContext2D;vtkViewsInfovis
endif(VIENNALS_STATIC_BUILD)
list(APPEND VIENNALS_LIBRARIES ${VTK_LIBRARIES})
list(APPEND VIENNALS_PYTHON_LIBRARIES ${VTK_LIBRARIES})
Expand Down
89 changes: 48 additions & 41 deletions Examples/AirGapDeposition/AirGapDeposition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,22 +16,26 @@
then grown directionally on top. \example AirGapDeposition.cpp
*/

using NumericType = float;

// implement own velocity field
class velocityField : public lsVelocityField<double> {
class velocityField : public lsVelocityField<NumericType> {
public:
double getScalarVelocity(const std::array<double, 3> & /*coordinate*/,
int /*material*/,
const std::array<double, 3> &normalVector) {
NumericType
getScalarVelocity(const std::array<NumericType, 3> & /*coordinate*/,
int /*material*/,
const std::array<NumericType, 3> &normalVector) {
// velocity is proportional to the normal vector
double velocity = std::abs(normalVector[0]) + std::abs(normalVector[1]);
NumericType velocity =
std::abs(normalVector[0]) + std::abs(normalVector[1]);
return velocity;
}

std::array<double, 3>
getVectorVelocity(const std::array<double, 3> & /*coordinate*/,
std::array<NumericType, 3>
getVectorVelocity(const std::array<NumericType, 3> & /*coordinate*/,
int /*material*/,
const std::array<double, 3> & /*normalVector*/) {
return std::array<double, 3>({});
const std::array<NumericType, 3> & /*normalVector*/) {
return std::array<NumericType, 3>({});
}
};

Expand All @@ -40,53 +44,55 @@ int main() {
constexpr int D = 2;
omp_set_num_threads(2);

double extent = 30;
double gridDelta = 0.5;
NumericType extent = 30;
NumericType gridDelta = 0.5;

double bounds[2 * D] = {-extent, extent, -extent, extent};
lsDomain<double, D>::BoundaryType boundaryCons[D];
boundaryCons[0] = lsDomain<double, D>::BoundaryType::REFLECTIVE_BOUNDARY;
boundaryCons[1] = lsDomain<double, D>::BoundaryType::INFINITE_BOUNDARY;
hrleCoordType bounds[2 * D] = {-extent, extent, -extent, extent};
lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
boundaryCons[0] = lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
boundaryCons[1] = lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;

auto substrate =
lsSmartPointer<lsDomain<double, D>>::New(bounds, boundaryCons, gridDelta);
auto substrate = lsSmartPointer<lsDomain<NumericType, D>>::New(
bounds, boundaryCons, gridDelta);

double origin[2] = {0., 0.};
double planeNormal[2] = {0., 1.};
NumericType origin[2] = {0., 0.};
NumericType planeNormal[2] = {0., 1.};

{
auto plane = lsSmartPointer<lsPlane<double, D>>::New(origin, planeNormal);
lsMakeGeometry<double, D>(substrate, plane).apply();
auto plane =
lsSmartPointer<lsPlane<NumericType, D>>::New(origin, planeNormal);
lsMakeGeometry<NumericType, D>(substrate, plane).apply();
}

{
std::cout << "Extracting..." << std::endl;
auto mesh = lsSmartPointer<lsMesh>::New();
lsToSurfaceMesh<double, D>(substrate, mesh).apply();
lsVTKWriter(mesh, "plane.vtk").apply();
auto mesh = lsSmartPointer<lsMesh<NumericType>>::New();
lsToSurfaceMesh<NumericType, D>(substrate, mesh).apply();
lsVTKWriter<NumericType>(mesh, "plane.vtk").apply();
}

{
// create layer used for booling
std::cout << "Creating box..." << std::endl;
auto trench = lsSmartPointer<lsDomain<double, D>>::New(bounds, boundaryCons,
gridDelta);
double minCorner[D] = {-extent / 6., -25.};
double maxCorner[D] = {extent / 6., 1.};
auto box = lsSmartPointer<lsBox<double, D>>::New(minCorner, maxCorner);
lsMakeGeometry<double, D>(trench, box).apply();
auto trench = lsSmartPointer<lsDomain<NumericType, D>>::New(
bounds, boundaryCons, gridDelta);
NumericType xlimit = extent / 6.;
NumericType minCorner[D] = {-xlimit, -25.};
NumericType maxCorner[D] = {xlimit, 1.};
auto box = lsSmartPointer<lsBox<NumericType, D>>::New(minCorner, maxCorner);
lsMakeGeometry<NumericType, D>(trench, box).apply();

{
std::cout << "Extracting..." << std::endl;
auto mesh = lsSmartPointer<lsMesh>::New();
lsToMesh<double, D>(trench, mesh).apply();
lsVTKWriter(mesh, "box.vtk").apply();
auto mesh = lsSmartPointer<lsMesh<NumericType>>::New();
lsToMesh<NumericType, D>(trench, mesh).apply();
lsVTKWriter<NumericType>(mesh, "box.vtk").apply();
}

// Create trench geometry
std::cout << "Booling trench..." << std::endl;
lsBooleanOperation<double, D>(substrate, trench,
lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
lsBooleanOperation<NumericType, D>(
substrate, trench, lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
.apply();
}

Expand All @@ -95,12 +101,12 @@ int main() {
// create new levelset for new material, which will be grown
// since it has to wrap around the substrate, just copy it
std::cout << "Creating new layer..." << std::endl;
auto newLayer = lsSmartPointer<lsDomain<double, D>>::New(substrate);
auto newLayer = lsSmartPointer<lsDomain<NumericType, D>>::New(substrate);

auto velocities = lsSmartPointer<velocityField>::New();

std::cout << "Advecting" << std::endl;
lsAdvect<double, D> advectionKernel;
lsAdvect<NumericType, D> advectionKernel;

// the level set to be advected has to be inserted last
// the other could be taken as a mask layer for advection
Expand All @@ -113,17 +119,18 @@ int main() {
// Now advect the level set 50 times, outputting every
// advection step. Save the physical time that
// passed during the advection.
double passedTime = 0.;
NumericType passedTime = 0.;
unsigned numberOfSteps = 60;
for (unsigned i = 0; i < numberOfSteps; ++i) {
advectionKernel.apply();
passedTime += advectionKernel.getAdvectedTime();

std::cout << "\rAdvection step " + std::to_string(i) + " / "
<< numberOfSteps << std::flush;
auto mesh = lsSmartPointer<lsMesh>::New();
lsToSurfaceMesh<double, D>(newLayer, mesh).apply();
lsVTKWriter(mesh, "trench" + std::to_string(i) + ".vtk").apply();
auto mesh = lsSmartPointer<lsMesh<NumericType>>::New();
lsToSurfaceMesh<NumericType, D>(newLayer, mesh).apply();
lsVTKWriter<NumericType>(mesh, "trench" + std::to_string(i) + ".vtk")
.apply();
}
std::cout << std::endl;
std::cout << "Time passed during advection: " << passedTime << std::endl;
Expand Down
94 changes: 50 additions & 44 deletions Examples/Deposition/Deposition.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,26 +17,26 @@
\example Deposition.cpp
*/

using NumericType = float;

// implement own velocity field
class velocityField : public lsVelocityField<double> {
class velocityField : public lsVelocityField<NumericType> {
public:
double
getScalarVelocity(const std::array<double, 3> & /*coordinate*/,
int /*material*/,
const std::array<double, 3>
& /*normalVector = hrleVectorType<double, 3>(0.)*/) {
NumericType getScalarVelocity(
const std::array<NumericType, 3> & /*coordinate*/, int /*material*/,
const std::array<NumericType, 3>
& /*normalVector = hrleVectorType<NumericType, 3>(0.)*/) {
// Some arbitrary velocity function of your liking
// (try changing it and see what happens :)
double velocity = 1.;
NumericType velocity = 1.;
return velocity;
}

std::array<double, 3>
getVectorVelocity(const std::array<double, 3> & /*coordinate*/,
int /*material*/,
const std::array<double, 3>
& /*normalVector = hrleVectorType<double, 3>(0.)*/) {
return std::array<double, 3>({}); // initialise to zero
std::array<NumericType, 3> getVectorVelocity(
const std::array<NumericType, 3> & /*coordinate*/, int /*material*/,
const std::array<NumericType, 3>
& /*normalVector = hrleVectorType<NumericType, 3>(0.)*/) {
return std::array<NumericType, 3>({}); // initialise to zero
}
};

Expand All @@ -45,58 +45,61 @@ int main() {
constexpr int D = 3;
omp_set_num_threads(4);

double extent = 30;
double gridDelta = 0.5;
NumericType extent = 30;
NumericType gridDelta = 0.5;

double bounds[2 * D] = {-extent, extent, -extent, extent, -extent, extent};
lsDomain<double, D>::BoundaryType boundaryCons[D];
lsDomain<NumericType, D>::BoundaryType boundaryCons[D];
for (unsigned i = 0; i < D - 1; ++i)
boundaryCons[i] = lsDomain<double, D>::BoundaryType::REFLECTIVE_BOUNDARY;
boundaryCons[2] = lsDomain<double, D>::BoundaryType::INFINITE_BOUNDARY;
boundaryCons[i] =
lsDomain<NumericType, D>::BoundaryType::REFLECTIVE_BOUNDARY;
boundaryCons[2] = lsDomain<NumericType, D>::BoundaryType::INFINITE_BOUNDARY;

auto substrate =
lsSmartPointer<lsDomain<double, D>>::New(bounds, boundaryCons, gridDelta);
auto substrate = lsSmartPointer<lsDomain<NumericType, D>>::New(
bounds, boundaryCons, gridDelta);

double origin[3] = {0., 0., 0.};
double planeNormal[3] = {0., 0., 1.};
NumericType origin[3] = {0., 0., 0.};
NumericType planeNormal[3] = {0., 0., 1.};

{
auto plane = lsSmartPointer<lsPlane<double, D>>::New(origin, planeNormal);
lsMakeGeometry<double, D>(substrate, plane).apply();
auto plane =
lsSmartPointer<lsPlane<NumericType, D>>::New(origin, planeNormal);
lsMakeGeometry<NumericType, D>(substrate, plane).apply();
}

{
auto trench = lsSmartPointer<lsDomain<double, D>>::New(bounds, boundaryCons,
gridDelta);
auto trench = lsSmartPointer<lsDomain<NumericType, D>>::New(
bounds, boundaryCons, gridDelta);
// make -x and +x greater than domain for numerical stability
double minCorner[D] = {-extent - 1, -extent / 4., -15.};
double maxCorner[D] = {extent + 1, extent / 4., 1.};
auto box = lsSmartPointer<lsBox<double, D>>::New(minCorner, maxCorner);
lsMakeGeometry<double, D>(trench, box).apply();
NumericType ylimit = extent / 4.;
NumericType minCorner[D] = {-extent - 1, -ylimit, -15.};
NumericType maxCorner[D] = {extent + 1, ylimit, 1.};
auto box = lsSmartPointer<lsBox<NumericType, D>>::New(minCorner, maxCorner);
lsMakeGeometry<NumericType, D>(trench, box).apply();

// Create trench geometry
lsBooleanOperation<double, D>(substrate, trench,
lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
lsBooleanOperation<NumericType, D>(
substrate, trench, lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
.apply();
}

{
std::cout << "Extracting..." << std::endl;
auto mesh = lsSmartPointer<lsMesh>::New();
lsToSurfaceMesh<double, D>(substrate, mesh).apply();
lsVTKWriter(mesh, "trench-0.vtk").apply();
auto mesh = lsSmartPointer<lsMesh<NumericType>>::New();
lsToSurfaceMesh<NumericType, D>(substrate, mesh).apply();
lsVTKWriter<NumericType>(mesh, "trench-0.vtk").apply();
}

// Now grow new material isotropically

// create new levelset for new material, which will be grown
// since it has to wrap around the substrate, just copy it
auto newLayer = lsSmartPointer<lsDomain<double, D>>::New(substrate);
auto newLayer = lsSmartPointer<lsDomain<NumericType, D>>::New(substrate);

auto velocities = lsSmartPointer<velocityField>::New();

std::cout << "Advecting" << std::endl;
lsAdvect<double, D> advectionKernel;
lsAdvect<NumericType, D> advectionKernel;

// the level set to be advected has to be inserted last
// the other could be taken as a mask layer for advection
Expand All @@ -106,20 +109,23 @@ int main() {
advectionKernel.setVelocityField(velocities);
// advectionKernel.setAdvectionTime(4.);
unsigned counter = 1;
for (double time = 0; time < 4.; time += advectionKernel.getAdvectedTime()) {
for (NumericType time = 0; time < 4.;
time += advectionKernel.getAdvectedTime()) {
advectionKernel.apply();

auto mesh = lsSmartPointer<lsMesh>::New();
lsToSurfaceMesh<double, D>(newLayer, mesh).apply();
lsVTKWriter(mesh, "trench-" + std::to_string(counter) + ".vtk").apply();
auto mesh = lsSmartPointer<lsMesh<NumericType>>::New();
lsToSurfaceMesh<NumericType, D>(newLayer, mesh).apply();
lsVTKWriter<NumericType>(mesh, "trench-" + std::to_string(counter) + ".vtk")
.apply();

lsToMesh<double, D>(newLayer, mesh).apply();
lsVTKWriter(mesh, "LS-" + std::to_string(counter) + ".vtk").apply();
lsToMesh<NumericType, D>(newLayer, mesh).apply();
lsVTKWriter<NumericType>(mesh, "LS-" + std::to_string(counter) + ".vtk")
.apply();

++counter;
}

// double advectionSteps = advectionKernel.getNumberOfTimeSteps();
// NumericType advectionSteps = advectionKernel.getNumberOfTimeSteps();
// std::cout << "Number of Advection steps taken: " << advectionSteps
// << std::endl;

Expand Down
Loading

0 comments on commit f58b429

Please sign in to comment.