zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry
2012-04-04 17:44:47 +01:00
parent fca4d338ea e8b03e3a25
commit e0d7783f00
24 changed files with 1155 additions and 7671 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/test/tetTetOverlap/Make/files Normal file
View File

@ -0,0 +1,3 @@
Test-tetTetOverlap.C
EXE = $(FOAM_USER_APPBIN)/Test-tetTetOverlap

5
applications/test/tetTetOverlap/Make/options Normal file
View File

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lmeshTools

164
applications/test/tetTetOverlap/Test-tetTetOverlap.C Normal file
View File

@ -0,0 +1,164 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
Test-tetTetOverlap
Description
Overlap volume of two tets
\*---------------------------------------------------------------------------*/
#include "tetrahedron.H"
#include "OFstream.H"
#include "meshTools.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void writeOBJ
(
Ostream& os,
label& vertI,
const tetPoints& tet
)
{
forAll(tet, fp)
{
meshTools::writeOBJ(os, tet[fp]);
}
os << "l " << vertI+1 << ' ' << vertI+2 << nl
<< "l " << vertI+1 << ' ' << vertI+3 << nl
<< "l " << vertI+1 << ' ' << vertI+4 << nl
<< "l " << vertI+2 << ' ' << vertI+3 << nl
<< "l " << vertI+2 << ' ' << vertI+4 << nl
<< "l " << vertI+3 << ' ' << vertI+4 << nl;
vertI += 4;
}
int main(int argc, char *argv[])
{
tetPoints A
(
point(0, 0, 0),
point(1, 0, 0),
point(1, 1, 0),
point(1, 1, 1)
);
const tetPointRef tetA = A.tet();
tetPoints B
(
point(0.1, 0.1, 0.1),
point(1.1, 0.1, 0.1),
point(1.1, 1.1, 0.1),
point(1.1, 1.1, 1.1)
);
const tetPointRef tetB = B.tet();
tetPointRef::tetIntersectionList insideTets;
label nInside = 0;
tetPointRef::tetIntersectionList outsideTets;
label nOutside = 0;
tetA.tetOverlap
(
tetB,
insideTets,
nInside,
outsideTets,
nOutside
);
// Dump to file
// ~~~~~~~~~~~~
{
OFstream str("tetA.obj");
Info<< "Writing A to " << str.name() << endl;
label vertI = 0;
writeOBJ(str, vertI, A);
}
{
OFstream str("tetB.obj");
Info<< "Writing B to " << str.name() << endl;
label vertI = 0;
writeOBJ(str, vertI, B);
}
{
OFstream str("inside.obj");
Info<< "Writing parts of A inside B to " << str.name() << endl;
label vertI = 0;
for (label i = 0; i < nInside; ++i)
{
writeOBJ(str, vertI, insideTets[i]);
}
}
{
OFstream str("outside.obj");
Info<< "Writing parts of A outside B to " << str.name() << endl;
label vertI = 0;
for (label i = 0; i < nOutside; ++i)
{
writeOBJ(str, vertI, outsideTets[i]);
}
}
// Check
// ~~~~~
Info<< "Vol A:" << tetA.mag() << endl;
scalar volInside = 0;
for (label i = 0; i < nInside; ++i)
{
volInside += insideTets[i].tet().mag();
}
Info<< "Vol A inside B:" << volInside << endl;
scalar volOutside = 0;
for (label i = 0; i < nOutside; ++i)
{
volOutside += outsideTets[i].tet().mag();
}
Info<< "Vol A outside B:" << volOutside << endl;
Info<< "Sum inside and outside:" << volInside+volOutside << endl;
if (mag(volInside+volOutside-tetA.mag()) > SMALL)
{
FatalErrorIn("Test-tetetOverlap")
<< "Tet volumes do not sum up to input tet."
<< exit(FatalError);
}
return 0;
}
// ************************************************************************* //

10
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Allwmake Executable file
View File

@ -0,0 +1,10 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
if [ -d "${FASTDUALOCTREE_SRC_PATH}" ]
then
wmake
fi
# ----------------------------------------------------------------- end-of-file

20
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/files
View File

@ -1,6 +1,22 @@
cvMeshSurfaceSimplify.C /*
cvMeshSurfaceSimplify_non_octree.C
MarchingCubes/MarchingCubes.cpp MarchingCubes/MarchingCubes.cpp
MarchingCubes/ply.c MarchingCubes/ply.c
*/
/*
MarchingCubes = fastdualoctree_sgp
$(MarchingCubes)/data_access.cpp
$(MarchingCubes)/fparser.cpp
$(MarchingCubes)/fpoptimizer.cpp
$(MarchingCubes)/MarchingCubes.cpp
$(MarchingCubes)/mc_draw.cpp
$(MarchingCubes)/morton.cpp
$(MarchingCubes)/opt_octree.cpp
$(MarchingCubes)/hash_octree.cpp
*/
cvMeshSurfaceSimplify.C
EXE = $(FOAM_APPBIN)/cvMeshSurfaceSimplify EXE = $(FOAM_APPBIN)/cvMeshSurfaceSimplify

9
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/Make/options
View File

@ -1,7 +1,12 @@
MarchingCubes = fastdualoctree_sgp
include $(GENERAL_RULES)/CGAL include $(GENERAL_RULES)/CGAL
EXE_INC = \ EXE_INC = \
-IMarchingCubes \ -DUNIX \
-Wno-old-style-cast \
/* -IMarchingCubes */ \
-I$(FASTDUALOCTREE_SRC_PATH) \
-I../conformalVoronoiMesh/lnInclude \ -I../conformalVoronoiMesh/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \ -I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude \ -I$(LIB_SRC)/triSurface/lnInclude \
@ -9,6 +14,8 @@ EXE_INC = \
EXE_LIBS = \ EXE_LIBS = \
$(CGAL_LIBS) \ $(CGAL_LIBS) \
-L$(FASTDUALOCTREE_SRC_PATH) -lperf_main \
-lGL \
-lconformalVoronoiMesh \ -lconformalVoronoiMesh \
-ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \ -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lscotchDecomp \
-ledgeMesh \ -ledgeMesh \

2322
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/LookUpTable.h
View File

File diff suppressed because it is too large Load Diff

1302
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/MarchingCubes.cpp
View File

File diff suppressed because it is too large Load Diff

343
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/MarchingCubes.h
View File

@ -1,343 +0,0 @@
/**
* @file MarchingCubes.h
* @author Thomas Lewiner <thomas.lewiner@polytechnique.org>
* @author Math Dept, PUC-Rio
* @version 0.2
* @date 12/08/2002
*
* @brief MarchingCubes Algorithm
*/
//________________________________________________
#ifndef _MARCHINGCUBES_H_
#define _MARCHINGCUBES_H_
#if !defined(WIN32) || defined(__CYGWIN__)
#pragma interface
#endif // WIN32
//_____________________________________________________________________________
// types
/** unsigned char alias */
typedef unsigned char uchar ;
/** signed char alias */
typedef signed char schar ;
/** isovalue alias */
typedef float real ;
//-----------------------------------------------------------------------------
// Vertex structure
/** \struct Vertex "MarchingCubes.h" MarchingCubes
* Position and normal of a vertex
* \brief vertex structure
* \param x X coordinate
* \param y Y coordinate
* \param z Z coordinate
* \param nx X component of the normal
* \param ny Y component of the normal
* \param nz Z component of the normal
*/
typedef struct
{
real x, y, z ; /**< Vertex coordinates */
real nx, ny, nz ; /**< Vertex normal */
} Vertex ;
//-----------------------------------------------------------------------------
// Triangle structure
/** \struct Triangle "MarchingCubes.h" MarchingCubes
* Indices of the oriented triange vertices
* \brief triangle structure
* \param v1 First vertex index
* \param v2 Second vertex index
* \param v3 Third vertex index
*/
typedef struct
{
int v1,v2,v3 ; /**< Triangle vertices */
} Triangle ;
//_____________________________________________________________________________
//_____________________________________________________________________________
/** Marching Cubes algorithm wrapper */
/** \class MarchingCubes
* \brief Marching Cubes algorithm.
*/
class MarchingCubes
//-----------------------------------------------------------------------------
{
// Constructors
public :
/**
* Main and default constructor
* \brief constructor
* \param size_x width of the grid
* \param size_y depth of the grid
* \param size_z height of the grid
*/
MarchingCubes ( const int size_x = -1, const int size_y = -1, const int size_z = -1 ) ;
/** Destructor */
~MarchingCubes() ;
//-----------------------------------------------------------------------------
// Accessors
public :
/** accesses the number of vertices of the generated mesh */
inline const int nverts() const { return _nverts ; }
/** accesses the number of triangles of the generated mesh */
inline const int ntrigs() const { return _ntrigs ; }
/** accesses a specific vertex of the generated mesh */
inline Vertex * vert( const int i ) const { if( i < 0 || i >= _nverts ) return ( Vertex *)NULL ; return _vertices + i ; }
/** accesses a specific triangle of the generated mesh */
inline Triangle * trig( const int i ) const { if( i < 0 || i >= _ntrigs ) return (Triangle*)NULL ; return _triangles + i ; }
/** accesses the vertex buffer of the generated mesh */
inline Vertex *vertices () { return _vertices ; }
/** accesses the triangle buffer of the generated mesh */
inline Triangle *triangles() { return _triangles ; }
/** accesses the width of the grid */
inline const int size_x() const { return _size_x ; }
/** accesses the depth of the grid */
inline const int size_y() const { return _size_y ; }
/** accesses the height of the grid */
inline const int size_z() const { return _size_z ; }
/**
* changes the size of the grid
* \param size_x width of the grid
* \param size_y depth of the grid
* \param size_z height of the grid
*/
inline void set_resolution( const int size_x, const int size_y, const int size_z ) { _size_x = size_x ; _size_y = size_y ; _size_z = size_z ; }
/**
* selects wether the algorithm will use the enhanced topologically controlled lookup table or the original MarchingCubes
* \param originalMC true for the original Marching Cubes
*/
inline void set_method ( const bool originalMC = false ) { _originalMC = originalMC ; }
/**
* selects to use data from another class
* \param data is the pointer to the external data, allocated as a size_x*size_y*size_z vector running in x first
*/
inline void set_ext_data ( real *data )
{ if( !_ext_data ) delete [] _data ; _ext_data = data != NULL ; if( _ext_data ) _data = data ; }
/**
* selects to allocate data
*/
inline void set_int_data () { _ext_data = false ; _data = NULL ; }
// Data access
/**
* accesses a specific cube of the grid
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline const real get_data ( const int i, const int j, const int k ) const { return _data[ i + j*_size_x + k*_size_x*_size_y] ; }
/**
* sets a specific cube of the grid
* \param val new value for the cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline void set_data ( const real val, const int i, const int j, const int k ) { _data[ i + j*_size_x + k*_size_x*_size_y] = val ; }
// Data initialization
/** inits temporary structures (must set sizes before call) : the grid and the vertex index per cube */
void init_temps () ;
/** inits all structures (must set sizes before call) : the temporary structures and the mesh buffers */
void init_all () ;
/** clears temporary structures : the grid and the main */
void clean_temps() ;
/** clears all structures : the temporary structures and the mesh buffers */
void clean_all () ;
//-----------------------------------------------------------------------------
// Exportation
public :
/**
* PLY exportation of the generated mesh
* \param fn name of the PLY file to create
* \param bin if true, the PLY will be written in binary mode
*/
void writePLY( const char *fn, bool bin = false ) ;
/**
* PLY importation of a mesh
* \param fn name of the PLY file to read from
*/
void readPLY( const char *fn ) ;
/**
* VRML / Open Inventor exportation of the generated mesh
* \param fn name of the IV file to create
*/
void writeIV ( const char *fn ) ;
/**
* ISO exportation of the input grid
* \param fn name of the ISO file to create
*/
void writeISO( const char *fn ) ;
//-----------------------------------------------------------------------------
// Algorithm
public :
/**
* Main algorithm : must be called after init_all
* \param iso isovalue
*/
void run( real iso = (real)0.0 ) ;
protected :
/** tesselates one cube */
void process_cube () ;
/** tests if the components of the tesselation of the cube should be connected by the interior of an ambiguous face */
bool test_face ( schar face ) ;
/** tests if the components of the tesselation of the cube should be connected through the interior of the cube */
bool test_interior( schar s ) ;
//-----------------------------------------------------------------------------
// Operations
protected :
/**
* computes almost all the vertices of the mesh by interpolation along the cubes edges
* \param iso isovalue
*/
void compute_intersection_points( real iso ) ;
/**
* routine to add a triangle to the mesh
* \param trig the code for the triangle as a sequence of edges index
* \param n the number of triangles to produce
* \param v12 the index of the interior vertex to use, if necessary
*/
void add_triangle ( const char* trig, char n, int v12 = -1 ) ;
/** tests and eventually doubles the vertex buffer capacity for a new vertex insertion */
void test_vertex_addition() ;
/** adds a vertex on the current horizontal edge */
int add_x_vertex() ;
/** adds a vertex on the current longitudinal edge */
int add_y_vertex() ;
/** adds a vertex on the current vertical edge */
int add_z_vertex() ;
/** adds a vertex inside the current cube */
int add_c_vertex() ;
/**
* interpolates the horizontal gradient of the implicit function at the lower vertex of the specified cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
real get_x_grad( const int i, const int j, const int k ) const ;
/**
* interpolates the longitudinal gradient of the implicit function at the lower vertex of the specified cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
real get_y_grad( const int i, const int j, const int k ) const ;
/**
* interpolates the vertical gradient of the implicit function at the lower vertex of the specified cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
real get_z_grad( const int i, const int j, const int k ) const ;
/**
* accesses the pre-computed vertex index on the lower horizontal edge of a specific cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline int get_x_vert( const int i, const int j, const int k ) const { return _x_verts[ i + j*_size_x + k*_size_x*_size_y] ; }
/**
* accesses the pre-computed vertex index on the lower longitudinal edge of a specific cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline int get_y_vert( const int i, const int j, const int k ) const { return _y_verts[ i + j*_size_x + k*_size_x*_size_y] ; }
/**
* accesses the pre-computed vertex index on the lower vertical edge of a specific cube
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline int get_z_vert( const int i, const int j, const int k ) const { return _z_verts[ i + j*_size_x + k*_size_x*_size_y] ; }
/**
* sets the pre-computed vertex index on the lower horizontal edge of a specific cube
* \param val the index of the new vertex
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline void set_x_vert( const int val, const int i, const int j, const int k ) { _x_verts[ i + j*_size_x + k*_size_x*_size_y] = val ; }
/**
* sets the pre-computed vertex index on the lower longitudinal edge of a specific cube
* \param val the index of the new vertex
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline void set_y_vert( const int val, const int i, const int j, const int k ) { _y_verts[ i + j*_size_x + k*_size_x*_size_y] = val ; }
/**
* sets the pre-computed vertex index on the lower vertical edge of a specific cube
* \param val the index of the new vertex
* \param i abscisse of the cube
* \param j ordinate of the cube
* \param k height of the cube
*/
inline void set_z_vert( const int val, const int i, const int j, const int k ) { _z_verts[ i + j*_size_x + k*_size_x*_size_y] = val ; }
/** prints cube for debug */
void print_cube() ;
//-----------------------------------------------------------------------------
// Elements
protected :
bool _originalMC ; /**< selects wether the algorithm will use the enhanced topologically controlled lookup table or the original MarchingCubes */
bool _ext_data ; /**< selects wether to allocate data or use data from another class */
int _size_x ; /**< width of the grid */
int _size_y ; /**< depth of the grid */
int _size_z ; /**< height of the grid */
real *_data ; /**< implicit function values sampled on the grid */
int *_x_verts ; /**< pre-computed vertex indices on the lower horizontal edge of each cube */
int *_y_verts ; /**< pre-computed vertex indices on the lower longitudinal edge of each cube */
int *_z_verts ; /**< pre-computed vertex indices on the lower vertical edge of each cube */
int _nverts ; /**< number of allocated vertices in the vertex buffer */
int _ntrigs ; /**< number of allocated triangles in the triangle buffer */
int _Nverts ; /**< size of the vertex buffer */
int _Ntrigs ; /**< size of the triangle buffer */
Vertex *_vertices ; /**< vertex buffer */
Triangle *_triangles ; /**< triangle buffer */
int _i ; /**< abscisse of the active cube */
int _j ; /**< height of the active cube */
int _k ; /**< ordinate of the active cube */
real _cube[8] ; /**< values of the implicit function on the active cube */
uchar _lut_entry ; /**< cube sign representation in [0..255] */
uchar _case ; /**< case of the active cube in [0..15] */
uchar _config ; /**< configuration of the active cube */
uchar _subconfig ; /**< subconfiguration of the active cube */
};
//_____________________________________________________________________________
#endif // _MARCHINGCUBES_H_

3325
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/ply.c
View File

File diff suppressed because it is too large Load Diff

233
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/MarchingCubes/ply.h
View File

@ -1,233 +0,0 @@
/*
Header for PLY polygon files.
- Greg Turk
A PLY file contains a single polygonal _object_.
An object is composed of lists of _elements_. Typical elements are
vertices, faces, edges and materials.
Each type of element for a given object has one or more _properties_
associated with the element type. For instance, a vertex element may
have as properties three floating-point values x,y,z and three unsigned
chars for red, green and blue.
-----------------------------------------------------------------------
Copyright (c) 1998 Georgia Institute of Technology. All rights reserved.
Permission to use, copy, modify and distribute this software and its
documentation for any purpose is hereby granted without fee, provided
that the above copyright notice and this permission notice appear in
all copies of this software and that you do not sell the software.
THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
*/
#ifndef __PLY_H__
#define __PLY_H__
#ifdef __cplusplus
extern "C" {
#endif
#include <stdio.h>
#include <stddef.h>
#define PLY_ASCII 1 /* ascii PLY file */
#define PLY_BINARY_BE 2 /* binary PLY file, big endian */
#define PLY_BINARY_LE 3 /* binary PLY file, little endian */
#define PLY_OKAY 0 /* ply routine worked okay */
#define PLY_ERROR -1 /* error in ply routine */
/* scalar data types supported by PLY format */
#define StartType 0
#define Int8 1
#define Int16 2
#define Int32 3
#define Uint8 4
#define Uint16 5
#define Uint32 6
#define Float32 7
#define Float64 8
#define EndType 9
#define PLY_SCALAR 0
#define PLY_LIST 1
#define PLY_STRING 2
typedef struct PlyProperty { /* description of a property */
char *name; /* property name */
int external_type; /* file's data type */
int internal_type; /* program's data type */
int offset; /* offset bytes of prop in a struct */
int is_list; /* 0 = scalar, 1 = list, 2 = char string */
int count_external; /* file's count type */
int count_internal; /* program's count type */
int count_offset; /* offset byte for list count */
} PlyProperty;
typedef struct PlyElement { /* description of an element */
char *name; /* element name */
int num; /* number of elements in this object */
int size; /* size of element (bytes) or -1 if variable */
int nprops; /* number of properties for this element */
PlyProperty **props; /* list of properties in the file */
char *store_prop; /* flags: property wanted by user? */
int other_offset; /* offset to un-asked-for props, or -1 if none*/
int other_size; /* size of other_props structure */
} PlyElement;
typedef struct PlyOtherProp { /* describes other properties in an element */
char *name; /* element name */
int size; /* size of other_props */
int nprops; /* number of properties in other_props */
PlyProperty **props; /* list of properties in other_props */
} PlyOtherProp;
typedef struct OtherData { /* for storing other_props for an other element */
void *other_props;
} OtherData;
typedef struct OtherElem { /* data for one "other" element */
char *elem_name; /* names of other elements */
int elem_count; /* count of instances of each element */
OtherData **other_data; /* actual property data for the elements */
PlyOtherProp *other_props; /* description of the property data */
} OtherElem;
typedef struct PlyOtherElems { /* "other" elements, not interpreted by user */
int num_elems; /* number of other elements */
OtherElem *other_list; /* list of data for other elements */
} PlyOtherElems;
#define AVERAGE_RULE 1
#define MAJORITY_RULE 2
#define MINIMUM_RULE 3
#define MAXIMUM_RULE 4
#define SAME_RULE 5
#define RANDOM_RULE 6
typedef struct PlyPropRules { /* rules for combining "other" properties */
PlyElement *elem; /* element whose rules we are making */
int *rule_list; /* types of rules (AVERAGE_PLY, MAJORITY_PLY, etc.) */
int nprops; /* number of properties we're combining so far */
int max_props; /* maximum number of properties we have room for now */
void **props; /* list of properties we're combining */
float *weights; /* list of weights of the properties */
} PlyPropRules;
typedef struct PlyRuleList {
char *name; /* name of the rule */
char *element; /* name of element that rule applies to */
char *property; /* name of property that rule applies to */
struct PlyRuleList *next; /* pointer for linked list of rules */
} PlyRuleList;
typedef struct PlyFile { /* description of PLY file */
FILE *fp; /* file pointer */
int file_type; /* ascii or binary */
float version; /* version number of file */
int num_elem_types; /* number of element types of object */
PlyElement **elems; /* list of elements */
int num_comments; /* number of comments */
char **comments; /* list of comments */
int num_obj_info; /* number of items of object information */
char **obj_info; /* list of object info items */
PlyElement *which_elem; /* element we're currently reading or writing */
PlyOtherElems *other_elems; /* "other" elements from a PLY file */
PlyPropRules *current_rules; /* current propagation rules */
PlyRuleList *rule_list; /* rule list from user */
} PlyFile;
/* memory allocation */
/*
extern char *my_alloc();
*/
#define myalloc(mem_size) my_alloc((mem_size), __LINE__, __FILE__)
/* old routines */
#if 0
extern PlyFile *ply_write(FILE *, int, char **, int);
extern PlyFile *ply_read(FILE *, int *, char ***);
extern PlyFile *ply_open_for_reading( char *, int *, char ***, int *, float *);
extern void ply_close(PlyFile *);
extern PlyOtherProp *ply_get_other_properties(PlyFile *, char *, int);
#endif
extern void ply_describe_property( PlyFile * , char * , PlyProperty * );
extern void ply_get_property( PlyFile * , char * , PlyProperty * );
extern void ply_get_element( PlyFile * , void * );
/*** delcaration of routines ***/
PlyOtherElems *get_other_element_ply( PlyFile * );
PlyFile *read_ply( FILE * );
PlyFile *write_ply( FILE * , int, char ** , int );
extern PlyFile *open_for_writing_ply( char * , int, char ** , int );
void close_ply( PlyFile * );
void free_ply( PlyFile * );
void get_info_ply( PlyFile * , float * , int * );
void free_other_elements_ply( PlyOtherElems * );
void append_comment_ply( PlyFile * , char * );
void append_obj_info_ply( PlyFile * , char * );
void copy_comments_ply( PlyFile * , PlyFile * );
void copy_obj_info_ply( PlyFile * , PlyFile * );
char* *get_comments_ply( PlyFile * , int * );
char* *get_obj_info_ply( PlyFile * , int * );
char* *get_element_list_ply( PlyFile * , int * );
int setup_property_ply( PlyFile * , PlyProperty * );
void get_element_ply( PlyFile * , void * );
char *setup_element_read_ply( PlyFile * , int, int * );
PlyOtherProp *get_other_properties_ply( PlyFile * , int );
void element_count_ply( PlyFile * , char * , int );
void describe_element_ply( PlyFile * , char * , int );
void describe_property_ply( PlyFile * , PlyProperty * );
void describe_other_properties_ply( PlyFile * , PlyOtherProp * , int );
void describe_other_elements_ply( PlyFile * , PlyOtherElems * );
void get_element_setup_ply( PlyFile * , char * , int, PlyProperty * );
PlyProperty* *get_element_description_ply( PlyFile * , char * , int * , int * );
void element_layout_ply( PlyFile * , char * , int, int, PlyProperty * );
void header_complete_ply( PlyFile * );
void put_element_setup_ply( PlyFile * , char * );
void put_element_ply( PlyFile * , void * );
void put_other_elements_ply( PlyFile * );
PlyPropRules *init_rule_ply( PlyFile * , char * );
void modify_rule_ply( PlyPropRules * , char * , int );
void start_props_ply( PlyFile * , PlyPropRules * );
void weight_props_ply( PlyFile * , float, void * );
void *get_new_props_ply( PlyFile * );
void set_prop_rules_ply( PlyFile * , PlyRuleList * );
PlyRuleList *append_prop_rule( PlyRuleList * , char * , char * );
int matches_rule_name( char * );
int equal_strings( char * , char * );
char *recreate_command_line( int, char *argv[] );
#ifdef __cplusplus
}
#endif
#endif /* !__PLY_H__ */

497
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/cvMeshSurfaceSimplify.C
View File

@ -28,6 +28,7 @@ Description
Simplifies surfaces by resampling. Simplifies surfaces by resampling.
Uses Thomas Lewiner's topology preserving MarchingCubes. Uses Thomas Lewiner's topology preserving MarchingCubes.
(http://zeus.mat.puc-rio.br/tomlew/tomlew_uk.php)
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -37,13 +38,328 @@ Description
#include "conformationSurfaces.H" #include "conformationSurfaces.H"
#include "triSurfaceMesh.H" #include "triSurfaceMesh.H"
#include "MarchingCubes.h" #include "opt_octree.h"
#include "cube.h"
using namespace Foam; using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
class pointConversion
{
const vector scale_;
const vector offset_;
public:
//- Construct from components
pointConversion
(
const vector scale,
const vector offset
)
:
scale_(scale),
offset_(offset)
{}
inline Point toLocal(const Foam::point& pt) const
{
Foam::point p = cmptMultiply(scale_, (pt + offset_));
return Point(p.x(), p.y(), p.z());
}
inline Foam::point toGlobal(const Point& pt) const
{
point p(pt.x(), pt.y(), pt.z());
return cmptDivide(p, scale_) - offset_;
}
};
// For use in Fast-Dual Octree from Thomas Lewiner
class distanceCalc
:
public ::data_access
{
const Level min_level_;
const conformationSurfaces& geometryToConformTo_;
const pointConversion& converter_;
// Private Member Functions
scalar signedDistance(const Foam::point& pt) const
{
const searchableSurfaces& geometry = geometryToConformTo_.geometry();
const labelList& surfaces = geometryToConformTo_.surfaces();
static labelList nearestSurfaces;
static scalarField distance;
static pointField samples(1);
samples[0] = pt;
searchableSurfacesQueries::signedDistance
(
geometry,
surfaces,
samples,
scalarField(1, GREAT),
searchableSurface::OUTSIDE,
nearestSurfaces,
distance
);
return distance[0];
}
public:
// Constructors
//- Construct from components
distanceCalc
(
Level max_level_,
real iso_val_,
Level min_level,
const conformationSurfaces& geometryToConformTo,
const pointConversion& converter
)
:
data_access(max_level_,iso_val_),
min_level_(min_level),
geometryToConformTo_(geometryToConformTo),
converter_(converter)
{}
//- Destructor
virtual ~distanceCalc()
{}
// Member Functions
//- test function
virtual bool need_refine( const Cube &c )
{
int l = c.lv() ;
if( l >= _max_level ) return false;
if( l < min_level_ ) return true;
treeBoundBox bb
(
converter_.toGlobal
(
Point
(
c.xmin(),
c.ymin(),
c.zmin()
)
),
converter_.toGlobal
(
Point
(
c.xmax(),
c.ymax(),
c.zmax()
)
)
);
const searchableSurfaces& geometry =
geometryToConformTo_.geometry();
const labelList& surfaces =
geometryToConformTo_.surfaces();
//- Uniform refinement around surface
{
forAll(surfaces, i)
{
if (geometry[surfaces[i]].overlaps(bb))
{
return true;
}
}
return false;
}
////- Surface intersects bb (but not using intersection test)
//scalar ccDist = signedDistance(bb.midpoint());
//scalar ccVal = ccDist - _iso_val;
//if (mag(ccVal) < SMALL)
//{
// return true;
//}
//const pointField points(bb.points());
//forAll(points, pointI)
//{
// scalar pointVal = signedDistance(points[pointI]) - _iso_val;
// if (ccVal*pointVal < 0)
// {
// return true;
// }
//}
//return false;
////- Refinement based on intersection with multiple planes.
//// Does not work well - too high a ratio between
//// neighbouring cubes.
//const pointField points(bb.points());
//const edgeList& edges = treeBoundBox::edges;
//pointField start(edges.size());
//pointField end(edges.size());
//forAll(edges, i)
//{
// start[i] = points[edges[i][0]];
// end[i] = points[edges[i][1]];
//}
//Foam::List<Foam::List<pointIndexHit> > hitInfo;
//labelListList hitSurfaces;
//searchableSurfacesQueries::findAllIntersections
//(
// geometry,
// surfaces,
// start,
// end,
// hitSurfaces,
// hitInfo
//);
//
//// Count number of intersections
//label nInt = 0;
//forAll(hitSurfaces, edgeI)
//{
// nInt += hitSurfaces[edgeI].size();
//}
//
//if (nInt == 0)
//{
// // No surface intersected. See if there is one inside
// forAll(surfaces, i)
// {
// if (geometry[surfaces[i]].overlaps(bb))
// {
// return true;
// }
// }
// return false;
//}
//
//// Check multiple surfaces
//label baseSurfI = -1;
//forAll(hitSurfaces, edgeI)
//{
// const labelList& hSurfs = hitSurfaces[edgeI];
// forAll(hSurfs, i)
// {
// if (baseSurfI == -1)
// {
// baseSurfI = hSurfs[i];
// }
// else if (baseSurfI != hSurfs[i])
// {
// // Multiple surfaces
// return true;
// }
// }
//}
//
//// Get normals
//DynamicList<pointIndexHit> baseInfo(nInt);
//forAll(hitInfo, edgeI)
//{
// const Foam::List<pointIndexHit>& hits = hitInfo[edgeI];
// forAll(hits, i)
// {
// (void)hits[i].hitPoint();
// baseInfo.append(hits[i]);
// }
//}
//vectorField normals;
//geometry[surfaces[baseSurfI]].getNormal(baseInfo, normals);
//for (label i = 1; i < normals.size(); ++i)
//{
// if ((normals[0] & normals[i]) < 0.9)
// {
// return true;
// }
//}
//labelList regions;
//geometry[surfaces[baseSurfI]].getRegion(baseInfo, regions);
//for (label i = 1; i < regions.size(); ++i)
//{
// if (regions[0] != regions[i])
// {
// return true;
// }
//}
//return false;
//samples[0] = point(c.xmin(), c.ymin(), c.zmin());
//samples[1] = point(c.xmax(), c.ymin(), c.zmin());
//samples[2] = point(c.xmax(), c.ymax(), c.zmin());
//samples[3] = point(c.xmin(), c.ymax(), c.zmin());
//
//samples[4] = point(c.xmin(), c.ymin(), c.zmax());
//samples[5] = point(c.xmax(), c.ymin(), c.zmax());
//samples[6] = point(c.xmax(), c.ymax(), c.zmax());
//samples[7] = point(c.xmin(), c.ymax(), c.zmax());
//scalarField nearestDistSqr(8, GREAT);
//
//Foam::List<pointIndexHit> nearestInfo;
//surf_.findNearest(samples, nearestDistSqr, nearestInfo);
//vectorField normals;
//surf_.getNormal(nearestInfo, normals);
//
//for (label i = 1; i < normals.size(); ++i)
//{
// if ((normals[0] & normals[i]) < 0.5)
// {
// return true;
// }
//}
//return false;
//// Check if surface octree same level
//const labelList elems(surf_.tree().findBox(bb));
//
//if (elems.size() > 1)
//{
// return true;
//}
//else
//{
// return false;
//}
}
//- data function
virtual real value_at( const Cube &c )
{
return signedDistance(converter_.toGlobal(c)) - _iso_val;
}
};
// Main program: // Main program:
int main(int argc, char *argv[]) int main(int argc, char *argv[])
@ -52,19 +368,16 @@ int main(int argc, char *argv[])
( (
"Re-sample surfaces used in cvMesh operation" "Re-sample surfaces used in cvMesh operation"
); );
//argList::validArgs.append("inputFile");
argList::validArgs.append("(nx ny nz)");
argList::validArgs.append("outputName"); argList::validArgs.append("outputName");
#include "setRootCase.H" #include "setRootCase.H"
#include "createTime.H" #include "createTime.H"
runTime.functionObjects().off(); runTime.functionObjects().off();
const Vector<label> n(IStringStream(args.args()[1])()); const fileName exportName = args.args()[1];
const fileName exportName = args.args()[2];
Info<< "Reading surfaces as specified in the cvMeshDict and" Info<< "Reading surfaces as specified in the cvMeshDict and"
<< " writing re-sampled " << n << " to " << exportName << " writing a re-sampled surface to " << exportName
<< nl << endl; << nl << endl;
cpuTime timer; cpuTime timer;
@ -114,126 +427,100 @@ int main(int argc, char *argv[])
<< timer.cpuTimeIncrement() << " s." << nl << endl; << timer.cpuTimeIncrement() << " s." << nl << endl;
// Extend
treeBoundBox bb = geometryToConformTo.globalBounds();
{
const vector smallVec = 0.1*bb.span();
bb.min() -= smallVec;
bb.max() += smallVec;
}
Info<< "Meshing to bounding box " << bb << nl << endl;
const vector span(bb.span());
const vector d
(
span.x()/(n.x()-1),
span.y()/(n.y()-1),
span.z()/(n.z()-1)
);
MarchingCubes mc(span.x(), span.y(), span.z() ) ;
mc.set_resolution(n.x(), n.y(), n.z());
mc.init_all() ;
// Generate points
pointField points(mc.size_x()*mc.size_y()*mc.size_z());
label pointI = 0;
point pt;
for( int k = 0 ; k < mc.size_z() ; k++ )
{
pt.z() = bb.min().z() + k*d.z();
for( int j = 0 ; j < mc.size_y() ; j++ )
{
pt.y() = bb.min().y() + j*d.y();
for( int i = 0 ; i < mc.size_x() ; i++ )
{
pt.x() = bb.min().x() + i*d.x();
points[pointI++] = pt;
}
}
}
Info<< "Generated " << points.size() << " sampling points in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
// Compute data
const searchableSurfaces& geometry = geometryToConformTo.geometry(); const searchableSurfaces& geometry = geometryToConformTo.geometry();
const labelList& surfaces = geometryToConformTo.surfaces(); const labelList& surfaces = geometryToConformTo.surfaces();
scalarField signedDist;
labelList nearestSurfaces; const label minLevel = 2;
searchableSurfacesQueries::signedDistance
// The max cube size follows from the minLevel and the default cube size
// (1)
const scalar maxSize = 1.0 / (1 << minLevel);
const scalar halfMaxSize = 0.5*maxSize;
// Scale the geometry to fit within
// halfMaxSize .. 1-halfMaxSize
scalar wantedRange = 1.0-maxSize;
const treeBoundBox bb = geometryToConformTo.globalBounds();
const vector scale = cmptDivide
( (
geometry, point(wantedRange, wantedRange, wantedRange),
surfaces, bb.span()
points, );
scalarField(points.size(), sqr(GREAT)), const vector offset =
searchableSurface::OUTSIDE, // for non-closed surfaces treat as cmptDivide
// outside (
nearestSurfaces, point(halfMaxSize, halfMaxSize, halfMaxSize),
signedDist scale
)
-bb.min();
const pointConversion converter(scale, offset);
// Marching cubes
OptOctree octree;
distanceCalc ref
(
8, //maxLevel
0.0, //distance
minLevel, //minLevel
geometryToConformTo,
converter
); );
// Fill elements octree.refine(&ref);
pointI = 0; octree.set_impl(&ref);
for( int k = 0 ; k < mc.size_z() ; k++ )
{
for( int j = 0 ; j < mc.size_y() ; j++ )
{
for( int i = 0 ; i < mc.size_x() ; i++ )
{
mc.set_data(float(signedDist[pointI++]), i, j, k);
}
}
}
Info<< "Determined signed distance in = " Info<< "Calculated octree in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl; << timer.cpuTimeIncrement() << " s." << nl << endl;
MarchingCubes& mc = octree.mc();
mc.run() ; mc.clean_all() ;
octree.build_isosurface(&ref) ;
Info<< "Constructed iso surface in = " Info<< "Constructed iso surface of distance in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl; << timer.cpuTimeIncrement() << " s." << nl << endl;
mc.clean_temps() ;
// Write output file // Write output file
if (mc.ntrigs() > 0) if (mc.ntrigs() > 0)
{ {
Triangle* triangles = mc.triangles(); Triangle* triangles = mc.triangles();
List<labelledTri> tris(mc.ntrigs()); label nTris = mc.ntrigs();
forAll(tris, triI) Foam::DynamicList<labelledTri> tris(mc.ntrigs());
for (label triI = 0; triI < nTris; ++triI)
{ {
tris[triI] = labelledTri const Triangle& t = triangles[triI];
( if (t.v1 != t.v2 && t.v1 != t.v3 && t.v2 != t.v3)
triangles[triI].v1, {
triangles[triI].v2, tris.append
triangles[triI].v3, (
0 // region labelledTri
); (
triangles[triI].v1,
triangles[triI].v2,
triangles[triI].v3,
0 // region
)
);
}
} }
Vertex* vertices = mc.vertices(); Point* vertices = mc.vertices();
pointField points(mc.nverts()); pointField points(mc.nverts());
forAll(points, pointI) forAll(points, pointI)
{ {
Vertex& v = vertices[pointI]; const Point& v = vertices[pointI];
points[pointI] = point points[pointI] = converter.toGlobal(v);
(
bb.min().x() + v.x*span.x()/mc.size_x(),
bb.min().y() + v.y*span.y()/mc.size_y(),
bb.min().z() + v.z*span.z()/mc.size_z()
);
} }
@ -267,6 +554,7 @@ int main(int argc, char *argv[])
} }
triSurface s(tris, patches, points, true); triSurface s(tris, patches, points, true);
tris.clearStorage();
Info<< "Extracted triSurface in = " Info<< "Extracted triSurface in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl; << timer.cpuTimeIncrement() << " s." << nl << endl;
@ -274,7 +562,7 @@ int main(int argc, char *argv[])
// Find out region on local surface of nearest point // Find out region on local surface of nearest point
{ {
List<pointIndexHit> hitInfo; Foam::List<pointIndexHit> hitInfo;
labelList hitSurfaces; labelList hitSurfaces;
geometryToConformTo.findSurfaceNearest geometryToConformTo.findSurfaceNearest
( (
@ -342,7 +630,6 @@ int main(int argc, char *argv[])
mc.clean_all() ; mc.clean_all() ;
Info<< "End\n" << endl; Info<< "End\n" << endl;
return 0; return 0;

352
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/cvMeshSurfaceSimplify_non_octree.C Normal file
View File

@ -0,0 +1,352 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
cvMeshSurfaceSimplify
Description
Simplifies surfaces by resampling.
Uses Thomas Lewiner's topology preserving MarchingCubes.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "searchableSurfaces.H"
#include "conformationSurfaces.H"
#include "triSurfaceMesh.H"
#include "MarchingCubes.h"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
argList::addNote
(
"Re-sample surfaces used in cvMesh operation"
);
//argList::validArgs.append("inputFile");
argList::validArgs.append("(nx ny nz)");
argList::validArgs.append("outputName");
#include "setRootCase.H"
#include "createTime.H"
runTime.functionObjects().off();
const Vector<label> n(IStringStream(args.args()[1])());
const fileName exportName = args.args()[2];
Info<< "Reading surfaces as specified in the cvMeshDict and"
<< " writing re-sampled " << n << " to " << exportName
<< nl << endl;
cpuTime timer;
IOdictionary cvMeshDict
(
IOobject
(
"cvMeshDict",
runTime.system(),
runTime,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// Define/load all geometry
searchableSurfaces allGeometry
(
IOobject
(
"cvSearchableSurfaces",
runTime.constant(),
"triSurface",
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
cvMeshDict.subDict("geometry")
);
Info<< "Geometry read in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
Random rndGen(64293*Pstream::myProcNo());
conformationSurfaces geometryToConformTo
(
runTime,
rndGen,
allGeometry,
cvMeshDict.subDict("surfaceConformation")
);
Info<< "Set up geometry in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
// Extend
treeBoundBox bb = geometryToConformTo.globalBounds();
{
const vector smallVec = 0.1*bb.span();
bb.min() -= smallVec;
bb.max() += smallVec;
}
Info<< "Meshing to bounding box " << bb << nl << endl;
const vector span(bb.span());
const vector d
(
span.x()/(n.x()-1),
span.y()/(n.y()-1),
span.z()/(n.z()-1)
);
MarchingCubes mc(span.x(), span.y(), span.z() ) ;
mc.set_resolution(n.x(), n.y(), n.z());
mc.init_all() ;
// Generate points
pointField points(mc.size_x()*mc.size_y()*mc.size_z());
label pointI = 0;
point pt;
for( int k = 0 ; k < mc.size_z() ; k++ )
{
pt.z() = bb.min().z() + k*d.z();
for( int j = 0 ; j < mc.size_y() ; j++ )
{
pt.y() = bb.min().y() + j*d.y();
for( int i = 0 ; i < mc.size_x() ; i++ )
{
pt.x() = bb.min().x() + i*d.x();
points[pointI++] = pt;
}
}
}
Info<< "Generated " << points.size() << " sampling points in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
// Compute data
const searchableSurfaces& geometry = geometryToConformTo.geometry();
const labelList& surfaces = geometryToConformTo.surfaces();
scalarField signedDist;
labelList nearestSurfaces;
searchableSurfacesQueries::signedDistance
(
geometry,
surfaces,
points,
scalarField(points.size(), sqr(GREAT)),
searchableSurface::OUTSIDE, // for non-closed surfaces treat as
// outside
nearestSurfaces,
signedDist
);
// Fill elements
pointI = 0;
for( int k = 0 ; k < mc.size_z() ; k++ )
{
for( int j = 0 ; j < mc.size_y() ; j++ )
{
for( int i = 0 ; i < mc.size_x() ; i++ )
{
mc.set_data(float(signedDist[pointI++]), i, j, k);
}
}
}
Info<< "Determined signed distance in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
mc.run() ;
Info<< "Constructed iso surface in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
mc.clean_temps() ;
// Write output file
if (mc.ntrigs() > 0)
{
Triangle* triangles = mc.triangles();
List<labelledTri> tris(mc.ntrigs());
forAll(tris, triI)
{
tris[triI] = labelledTri
(
triangles[triI].v1,
triangles[triI].v2,
triangles[triI].v3,
0 // region
);
}
Vertex* vertices = mc.vertices();
pointField points(mc.nverts());
forAll(points, pointI)
{
Vertex& v = vertices[pointI];
points[pointI] = point
(
bb.min().x() + v.x*span.x()/mc.size_x(),
bb.min().y() + v.y*span.y()/mc.size_y(),
bb.min().z() + v.z*span.z()/mc.size_z()
);
}
// Find correspondence to original surfaces
labelList regionOffsets(surfaces.size());
label nRegions = 0;
forAll(surfaces, i)
{
const wordList& regions = geometry[surfaces[i]].regions();
regionOffsets[i] = nRegions;
nRegions += regions.size();
}
geometricSurfacePatchList patches(nRegions);
nRegions = 0;
forAll(surfaces, i)
{
const wordList& regions = geometry[surfaces[i]].regions();
forAll(regions, regionI)
{
patches[nRegions] = geometricSurfacePatch
(
"patch",
geometry[surfaces[i]].name() + "_" + regions[regionI],
nRegions
);
nRegions++;
}
}
triSurface s(tris, patches, points, true);
Info<< "Extracted triSurface in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
// Find out region on local surface of nearest point
{
List<pointIndexHit> hitInfo;
labelList hitSurfaces;
geometryToConformTo.findSurfaceNearest
(
s.faceCentres(),
scalarField(s.size(), sqr(GREAT)),
hitInfo,
hitSurfaces
);
// Get region
DynamicList<pointIndexHit> surfInfo(hitSurfaces.size());
DynamicList<label> surfIndices(hitSurfaces.size());
forAll(surfaces, surfI)
{
// Extract info on this surface
surfInfo.clear();
surfIndices.clear();
forAll(hitSurfaces, triI)
{
if (hitSurfaces[triI] == surfI)
{
surfInfo.append(hitInfo[triI]);
surfIndices.append(triI);
}
}
// Calculate sideness of these surface points
labelList region;
geometry[surfaces[surfI]].getRegion(surfInfo, region);
forAll(region, i)
{
label triI = surfIndices[i];
s[triI].region() = regionOffsets[surfI]+region[i];
}
}
}
Info<< "Re-patched surface in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
triSurfaceMesh smesh
(
IOobject
(
exportName,
runTime.constant(), // instance
"triSurface",
runTime, // registry
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
s
);
Info<< "writing surfMesh:\n "
<< smesh.searchableSurface::objectPath() << nl << endl;
smesh.searchableSurface::write();
Info<< "Written surface in = "
<< timer.cpuTimeIncrement() << " s." << nl << endl;
}
mc.clean_all() ;
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

BIN
applications/utilities/mesh/generation/cvMesh/cvMeshSurfaceSimplify/marching_cubes_jgt.zip
View File

Binary file not shown.

95
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -21,9 +21,6 @@ License
You should have received a copy of the GNU General Public License You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Calculation of shape function product for a tetrahedron
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "tetrahedron.H" #include "tetrahedron.H"
@ -32,6 +29,96 @@ Description
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Point, class PointRef>
void Foam::tetrahedron<Point, PointRef>::tetOverlap
(
const tetrahedron<Point, PointRef>& tetB,
tetIntersectionList& insideTets,
label& nInside,
tetIntersectionList& outsideTets,
label& nOutside
) const
{
// Work storage
tetIntersectionList cutInsideTets;
label nCutInside = 0;
nInside = 0;
storeOp inside(insideTets, nInside);
storeOp cutInside(cutInsideTets, nCutInside);
nOutside = 0;
storeOp outside(outsideTets, nOutside);
// Cut tetA with all inwards pointing faces of tetB. Any tets remaining
// in aboveTets are inside tetB.
{
// face0
plane pl0(tetB.b_, tetB.d_, tetB.c_);
// Cut and insert subtets into cutInsideTets (either by getting
// an index from freeSlots or by appending to insideTets) or
// insert into outsideTets
sliceWithPlane(pl0, cutInside, outside);
}
if (nCutInside == 0)
{
nInside = nCutInside;
return;
}
{
// face1
plane pl1(tetB.a_, tetB.c_, tetB.d_);
nInside = 0;
for (label i = 0; i < nCutInside; i++)
{
cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
}
if (nInside == 0)
{
return;
}
}
{
// face2
plane pl2(tetB.a_, tetB.d_, tetB.b_);
nCutInside = 0;
for (label i = 0; i < nInside; i++)
{
insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
}
if (nCutInside == 0)
{
nInside = nCutInside;
return;
}
}
{
// face3
plane pl3(tetB.a_, tetB.b_, tetB.c_);
nInside = 0;
for (label i = 0; i < nCutInside; i++)
{
cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
}
}
}
// (Probably very inefficient) minimum containment sphere calculation. // (Probably very inefficient) minimum containment sphere calculation.
// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf: // From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
// Sphere ctr is smallest one of // Sphere ctr is smallest one of

21
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H
View File

@ -87,6 +87,13 @@ class tetrahedron
{ {
public: public:
// Public typedefs
//- Storage type for tets originating from intersecting tets.
// (can possibly be smaller than 200)
typedef FixedList<tetPoints, 200> tetIntersectionList;
// Classes for use in sliceWithPlane. What to do with decomposition // Classes for use in sliceWithPlane. What to do with decomposition
// of tet. // of tet.
@ -111,11 +118,11 @@ public:
//- Store resulting tets //- Store resulting tets
class storeOp class storeOp
{ {
FixedList<tetPoints, 200>& tets_; tetIntersectionList& tets_;
label& nTets_; label& nTets_;
public: public:
inline storeOp(FixedList<tetPoints, 200>&, label&); inline storeOp(tetIntersectionList&, label&);
inline void operator()(const tetPoints&); inline void operator()(const tetPoints&);
}; };
@ -261,6 +268,16 @@ public:
BelowTetOp& belowOp BelowTetOp& belowOp
) const; ) const;
//- Decompose tet into tets inside and outside other tet
inline void tetOverlap
(
const tetrahedron<Point, PointRef>& tetB,
tetIntersectionList& insideTets,
label& nInside,
tetIntersectionList& outsideTets,
label& nOutside
) const;
//- Return (min)containment sphere, i.e. the smallest sphere with //- Return (min)containment sphere, i.e. the smallest sphere with
// all points inside. Returns pointHit with: // all points inside. Returns pointHit with:

2
src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H
View File

@ -556,7 +556,7 @@ inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
template<class Point, class PointRef> template<class Point, class PointRef>
inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
( (
FixedList<tetPoints, 200>& tets, tetIntersectionList& tets,
label& nTets label& nTets
) )
: :

32
src/finiteVolume/cfdTools/general/MRF/MRFZone.C
View File

@ -32,6 +32,7 @@ License
#include "faceSet.H" #include "faceSet.H"
#include "geometricOneField.H" #include "geometricOneField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::MRFZone, 0); defineTypeNameAndDebug(Foam::MRFZone, 0);
@ -273,8 +274,8 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
} }
} }
bool cellZoneFound = (cellZoneID_ != -1); bool cellZoneFound = (cellZoneID_ != -1);
reduce(cellZoneFound, orOp<bool>()); reduce(cellZoneFound, orOp<bool>());
if (!cellZoneFound) if (!cellZoneFound)
@ -307,6 +308,7 @@ void Foam::MRFZone::addCoriolis
const scalarField& V = mesh_.V(); const scalarField& V = mesh_.V();
vectorField& ddtUc = ddtU.internalField(); vectorField& ddtUc = ddtU.internalField();
const vectorField& Uc = U.internalField(); const vectorField& Uc = U.internalField();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
forAll(cells, i) forAll(cells, i)
@ -328,6 +330,7 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
const scalarField& V = mesh_.V(); const scalarField& V = mesh_.V();
vectorField& Usource = UEqn.source(); vectorField& Usource = UEqn.source();
const vectorField& U = UEqn.psi(); const vectorField& U = UEqn.psi();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
forAll(cells, i) forAll(cells, i)
@ -353,6 +356,7 @@ void Foam::MRFZone::addCoriolis
const scalarField& V = mesh_.V(); const scalarField& V = mesh_.V();
vectorField& Usource = UEqn.source(); vectorField& Usource = UEqn.source();
const vectorField& U = UEqn.psi(); const vectorField& U = UEqn.psi();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
forAll(cells, i) forAll(cells, i)
@ -368,6 +372,7 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const
const volVectorField& C = mesh_.C(); const volVectorField& C = mesh_.C();
const vector& origin = origin_.value(); const vector& origin = origin_.value();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
const labelList& cells = mesh_.cellZones()[cellZoneID_]; const labelList& cells = mesh_.cellZones()[cellZoneID_];
@ -395,7 +400,8 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const
{ {
label patchFacei = excludedFaces_[patchi][i]; label patchFacei = excludedFaces_[patchi][i];
U.boundaryField()[patchi][patchFacei] -= U.boundaryField()[patchi][patchFacei] -=
(Omega ^ (C.boundaryField()[patchi][patchFacei] - origin)); (Omega
^ (C.boundaryField()[patchi][patchFacei] - origin));
} }
} }
} }
@ -406,6 +412,7 @@ void Foam::MRFZone::absoluteVelocity(volVectorField& U) const
const volVectorField& C = mesh_.C(); const volVectorField& C = mesh_.C();
const vector& origin = origin_.value(); const vector& origin = origin_.value();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
const labelList& cells = mesh_.cellZones()[cellZoneID_]; const labelList& cells = mesh_.cellZones()[cellZoneID_];
@ -475,6 +482,7 @@ void Foam::MRFZone::absoluteFlux
void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
{ {
const vector& origin = origin_.value(); const vector& origin = origin_.value();
const vector& Omega = Omega_.value(); const vector& Omega = Omega_.value();
// Included patches // Included patches
@ -496,4 +504,24 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
} }
Foam::Ostream& Foam::operator<<(Ostream& os, const MRFZone& MRF)
{
os << indent << nl;
os.write(MRF.name_) << nl;
os << token::BEGIN_BLOCK << incrIndent << nl;
os.writeKeyword("origin") << MRF.origin_ << token::END_STATEMENT << nl;
os.writeKeyword("axis") << MRF.axis_ << token::END_STATEMENT << nl;
os.writeKeyword("omega") << MRF.omega_ << token::END_STATEMENT << nl;
if (MRF.excludedPatchNames_.size())
{
os.writeKeyword("nonRotatingPatches") << MRF.excludedPatchNames_
<< token::END_STATEMENT << nl;
}
os << decrIndent << token::END_BLOCK << nl;
return os;
}
// ************************************************************************* // // ************************************************************************* //

13
src/finiteVolume/cfdTools/general/MRF/MRFZone.H
View File

@ -48,6 +48,7 @@ SourceFiles
#include "surfaceFieldsFwd.H" #include "surfaceFieldsFwd.H"
#include "fvMatricesFwd.H" #include "fvMatricesFwd.H"
#include "fvMatrices.H" #include "fvMatrices.H"
#include "DataEntry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -87,7 +88,7 @@ class MRFZone
const dimensionedVector origin_; const dimensionedVector origin_;
dimensionedVector axis_; dimensionedVector axis_;
const dimensionedScalar omega_; dimensionedScalar omega_;
dimensionedVector Omega_; dimensionedVector Omega_;
@ -205,13 +206,11 @@ public:
void correctBoundaryVelocity(volVectorField& U) const; void correctBoundaryVelocity(volVectorField& U) const;
// Ostream Operator // IOstream operator
friend Ostream& operator<<(Ostream& os, const MRFZone& MRF);
friend Ostream& operator<<(Ostream& os, const MRFZone&)
{
notImplemented("Ostream& operator<<(Ostream& os, const MRFZone&)");
return os;
}
}; };

29
src/finiteVolume/cfdTools/general/MRF/MRFZones.C
View File

@ -49,8 +49,28 @@ Foam::MRFZones::MRFZones(const fvMesh& mesh)
IOobject::NO_WRITE IOobject::NO_WRITE
), ),
MRFZone::iNew(mesh) MRFZone::iNew(mesh)
),
mesh_(mesh)
{
if
(
Pstream::parRun()
&&
(
regIOobject::fileModificationChecking == timeStampMaster
|| regIOobject::fileModificationChecking == inotifyMaster
)
) )
{} {
WarningIn("MRFZones(const fvMesh&)")
<< "The MRFZones are not run time modifiable\n"
<< " using 'timeStampMaster' or 'inotifyMaster'\n"
<< " for the entry fileModificationChecking\n"
<< " in the etc/controlDict.\n"
<< " Use 'timeStamp' instead."
<< endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -161,4 +181,11 @@ void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
} }
bool Foam::MRFZones::readData(Istream& is)
{
PtrList<MRFZone>::read(is, MRFZone::iNew(mesh_));
return is.good();
}
// ************************************************************************* // // ************************************************************************* //

12
src/finiteVolume/cfdTools/general/MRF/MRFZones.H
View File

@ -52,6 +52,11 @@ class MRFZones
: :
public IOPtrList<MRFZone> public IOPtrList<MRFZone>
{ {
// Private data
//- Reference to mesh
const fvMesh& mesh_;
// Private Member Functions // Private Member Functions
@ -109,6 +114,13 @@ public:
//- Correct the boundary velocity for the roation of the MRF region //- Correct the boundary velocity for the roation of the MRF region
void correctBoundaryVelocity(volVectorField& U) const; void correctBoundaryVelocity(volVectorField& U) const;
// I-O
//- Read from Istream
virtual bool readData(Istream&);
}; };

19
src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -118,7 +118,7 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells()); inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells());
scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_; scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_;
labelListList& cellToCell = *cellToCellAddressingPtr_; const labelListList& cellToCell = cellToCellAddressing();
tetOverlapVolume overlapEngine; tetOverlapVolume overlapEngine;
@ -129,8 +129,8 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
if (overlapCells.size() > 0) if (overlapCells.size() > 0)
{ {
invVolCoeffs[celli].setSize(overlapCells.size()); invVolCoeffs[celli].setSize(overlapCells.size());
scalar v(0);
forAll (overlapCells, j) forAll(overlapCells, j)
{ {
label cellFrom = overlapCells[j]; label cellFrom = overlapCells[j];
treeBoundBox bbFromMesh treeBoundBox bbFromMesh
@ -142,7 +142,7 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
) )
); );
v = overlapEngine.cellCellOverlapVolumeMinDecomp scalar v = overlapEngine.cellCellOverlapVolumeMinDecomp
( (
toMesh_, toMesh_,
celli, celli,
@ -151,19 +151,14 @@ void Foam::meshToMesh::calculateInverseVolumeWeights() const
cellFrom, cellFrom,
bbFromMesh bbFromMesh
); );
invVolCoeffs[celli][j] = v/toMesh_.V()[celli]; invVolCoeffs[celli][j] = v/toMesh_.V()[celli];
}
if (celli == 2)
{
Info << "cellToCell :" << cellToCell[celli] << endl;
Info << "invVolCoeffs :" << invVolCoeffs[celli] << endl;
} }
} }
} }
} }
void Foam::meshToMesh::calculateCellToCellAddressing() const void Foam::meshToMesh::calculateCellToCellAddressing() const
{ {
if (debug) if (debug)
{ {

12
src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C
View File

@ -49,14 +49,14 @@ void Foam::tetOverlapVolume::tetTetOverlap
( (
const tetPoints& tetA, const tetPoints& tetA,
const tetPoints& tetB, const tetPoints& tetB,
FixedList<tetPoints, 200>& insideTets, tetIntersectionList& insideTets,
label& nInside, label& nInside,
FixedList<tetPoints, 200>& outsideTets, tetIntersectionList& outsideTets,
label& nOutside label& nOutside
) const ) const
{ {
// Work storage // Work storage
FixedList<tetPoints, 200> cutInsideTets; tetIntersectionList cutInsideTets;
label nCutInside = 0; label nCutInside = 0;
tetPointRef::storeOp inside(insideTets, nInside); tetPointRef::storeOp inside(insideTets, nInside);
@ -138,9 +138,9 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
const tetPoints& tetB const tetPoints& tetB
) const ) const
{ {
FixedList<tetPoints, 200> insideTets; tetIntersectionList insideTets;
label nInside = 0; label nInside = 0;
FixedList<tetPoints, 200> cutInsideTets; tetIntersectionList cutInsideTets;
label nCutInside = 0; label nCutInside = 0;
tetPointRef::storeOp inside(insideTets, nInside); tetPointRef::storeOp inside(insideTets, nInside);
@ -222,7 +222,7 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
const primitiveMesh& meshB, const primitiveMesh& meshB,
const label cellBI, const label cellBI,
const treeBoundBox& cellBbB const treeBoundBox& cellBbB
) ) const
{ {
const cell& cFacesA = meshA.cells()[cellAI]; const cell& cFacesA = meshA.cells()[cellAI];
const point& ccA = meshA.cellCentres()[cellAI]; const point& ccA = meshA.cellCentres()[cellAI];

6
src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H
View File

@ -60,9 +60,9 @@ class tetOverlapVolume
( (
const tetPoints& tetA, const tetPoints& tetA,
const tetPoints& tetB, const tetPoints& tetB,
FixedList<tetPoints, 200>& insideTets, tetIntersectionList& insideTets,
label& nInside, label& nInside,
FixedList<tetPoints, 200>& outsideTets, tetIntersectionList& outsideTets,
label& nOutside label& nOutside
) const; ) const;
@ -115,7 +115,7 @@ public:
const primitiveMesh& meshB, const primitiveMesh& meshB,
const label cellBI, const label cellBI,
const treeBoundBox& cellBbB const treeBoundBox& cellBbB
); ) const;
}; };