Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Kirill Terekhov
INMOST
Commits
cc18b29b
Commit
cc18b29b
authored
Feb 20, 2020
by
Kirill Terekhov
Browse files
Add various grid tools
parent
d3faa14f
Changes
15
Hide whitespace changes
Inline
Side-by-side
Examples/GridTools/CMakeLists.txt
View file @
cc18b29b
project
(
GridGen
)
add_executable
(
FixFaults fix_faults.cpp
)
add_executable
(
FixFaults test_fix_faults.cpp
)
add_executable
(
SetLayers test_set_layers.cpp
)
add_executable
(
MakeFaults test_make_faults.cpp
)
add_executable
(
FixTiny fix_tiny.cpp
)
add_executable
(
FixTinyCollapse fix_tiny_collapse.cpp
)
add_executable
(
UniteFaces unite_faces.cpp
)
add_executable
(
Dual dual.cpp
)
add_executable
(
Tetra tetra.cpp
)
add_executable
(
Slice slice.cpp
)
add_executable
(
SliceFunc slice_func
_test
.cpp
)
add_executable
(
SliceFunc
test_
slice_func.cpp
)
add_executable
(
SplitNonplanar split_nonplanar.cpp
)
add_executable
(
CollapseDegenerate collapse_degenerate.cpp
)
add_executable
(
Bnd2Stl bnd2stl.cpp
)
...
...
@@ -16,7 +18,7 @@ add_executable(Sew sew.cpp)
add_executable
(
Scale scale.cpp
)
add_executable
(
Move move.cpp
)
add_executable
(
Convert convert.cpp
)
add_executable
(
CubeTransform cube_transform.cpp
)
add_executable
(
CubeTransform
test_
cube_transform.cpp
)
add_executable
(
SameMeshDifference difference_same.cpp
)
add_executable
(
MatchSameMeshDifference difference_same_match.cpp
)
add_executable
(
MeshDifference difference_map.cpp
)
...
...
@@ -34,24 +36,61 @@ add_executable(sinusoidal_mesh sinusoidal_grid.cpp)
add_executable
(
split_faces split_faces.cpp
)
add_executable
(
glue_faces glue_faces.cpp
)
add_executable
(
check_collapse check_collapse.cpp
)
add_executable
(
test_f
racture test_fracture.cpp
)
add_executable
(
F
racture test_fracture.cpp
)
add_executable
(
segment_data segment_data.cpp
)
add_executable
(
ChangeTagName change_tag_name.cpp
)
add_library
(
FractureLib fracture.cpp fracture.h
)
add_library
(
SliceFuncLib slice_func.cpp slice_func.h
)
add_library
(
FixFaultsLib fix_faults.cpp fix_faults.h
)
add_library
(
CubeTransformLib cube_transform.cpp cube_transform.h
)
add_library
(
SetLayersLib set_layers.cpp set_layers.h
)
target_link_libraries
(
test_fracture inmost
)
target_link_libraries
(
test_fracture FractureLib
)
target_link_libraries
(
MakeFaults inmost
)
target_link_libraries
(
MakeFaults FractureLib
)
target_link_libraries
(
MakeFaults FixFaultsLib
)
target_link_libraries
(
MakeFaults SliceFuncLib
)
if
(
USE_MPI
)
message
(
"linking
test_fracture
with MPI"
)
target_link_libraries
(
test_fracture
${
MPI_LIBRARIES
}
)
message
(
"linking
MakeFaults
with MPI"
)
target_link_libraries
(
MakeFaults
${
MPI_LIBRARIES
}
)
if
(
MPI_LINK_FLAGS
)
set_target_properties
(
test_fracture
PROPERTIES LINK_FLAGS
"
${
MPI_LINK_FLAGS
}
"
)
set_target_properties
(
MakeFaults
PROPERTIES LINK_FLAGS
"
${
MPI_LINK_FLAGS
}
"
)
endif
()
endif
(
USE_MPI
)
install
(
TARGETS
test_fracture
EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
install
(
TARGETS
MakeFaults
EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
target_link_libraries
(
SetLayers inmost
)
target_link_libraries
(
SetLayers SetLayersLib
)
if
(
USE_MPI
)
message
(
"linking SetLayers with MPI"
)
target_link_libraries
(
SetLayers
${
MPI_LIBRARIES
}
)
if
(
MPI_LINK_FLAGS
)
set_target_properties
(
SetLayers PROPERTIES LINK_FLAGS
"
${
MPI_LINK_FLAGS
}
"
)
endif
()
endif
(
USE_MPI
)
install
(
TARGETS SetLayers EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
target_link_libraries
(
Fracture inmost
)
target_link_libraries
(
Fracture FractureLib
)
if
(
USE_MPI
)
message
(
"linking Fracture with MPI"
)
target_link_libraries
(
Fracture
${
MPI_LIBRARIES
}
)
if
(
MPI_LINK_FLAGS
)
set_target_properties
(
Fracture PROPERTIES LINK_FLAGS
"
${
MPI_LINK_FLAGS
}
"
)
endif
()
endif
(
USE_MPI
)
install
(
TARGETS Fracture EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
target_link_libraries
(
FixFaults inmost
)
target_link_libraries
(
FixFaults FixFaultsLib
)
if
(
USE_MPI
)
message
(
"linking FixFaults with MPI"
)
target_link_libraries
(
FixFaults
${
MPI_LIBRARIES
}
)
if
(
MPI_LINK_FLAGS
)
set_target_properties
(
FixFaults PROPERTIES LINK_FLAGS
"
${
MPI_LINK_FLAGS
}
"
)
endif
()
endif
(
USE_MPI
)
install
(
TARGETS FixFaults EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
target_link_libraries
(
FixTinyCollapse inmost
)
if
(
USE_MPI
)
...
...
@@ -95,16 +134,6 @@ endif(USE_MPI)
install
(
TARGETS glue_faces EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
target_link_libraries
(
FixFaults inmost
)
if
(
USE_MPI
)
message
(
"linking FixFaults with MPI"
)
target_link_libraries
(
FixFaults
${
MPI_LIBRARIES
}
)
if
(
MPI_LINK_FLAGS
)
set_target_properties
(
FixFaults PROPERTIES LINK_FLAGS
"
${
MPI_LINK_FLAGS
}
"
)
endif
()
endif
(
USE_MPI
)
install
(
TARGETS FixFaults EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools
)
target_link_libraries
(
FixTiny inmost
)
if
(
USE_MPI
)
message
(
"linking FixTiny with MPI"
)
...
...
@@ -267,6 +296,7 @@ install(TARGETS Convert EXPORT inmost-targets RUNTIME DESTINATION bin/GridTools)
target_link_libraries
(
CubeTransform inmost
)
target_link_libraries
(
CubeTransform CubeTransformLib
)
if
(
USE_MPI
)
message
(
"linking CubeTransform with MPI"
)
target_link_libraries
(
CubeTransform
${
MPI_LIBRARIES
}
)
...
...
@@ -444,8 +474,10 @@ install(TARGETS ChangeTagName EXPORT inmost-targets RUNTIME DESTINATION bin/Grid
set_property
(
TARGET FractureLib PROPERTY PUBLIC_HEADER
"
${
CMAKE_CURRENT_SOURCE_DIR
}
/fracture.h"
)
set_property
(
TARGET SliceFuncLib PROPERTY PUBLIC_HEADER
"
${
CMAKE_CURRENT_SOURCE_DIR
}
/slice_func.h"
)
set_property
(
TARGET FixFaultsLib PROPERTY PUBLIC_HEADER
"
${
CMAKE_CURRENT_SOURCE_DIR
}
/fix_faults.h"
)
set_property
(
TARGET CubeTransformLib PROPERTY PUBLIC_HEADER
"
${
CMAKE_CURRENT_SOURCE_DIR
}
/cube_transform.h"
)
install
(
TARGETS FractureLib SliceFuncLib EXPORT inmost-targets
install
(
TARGETS FractureLib SliceFuncLib
FixFaultsLib CubeTransformLib
EXPORT inmost-targets
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include
)
...
...
Examples/GridTools/cube_transform.cpp
View file @
cc18b29b
#include "inmost.h"
#include "cube_transform.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
std
::
string
mesh_name
=
"_TRANSFORMED"
;
using
namespace
INMOST
;
void
CubeTransform
::
PrintMapping
()
{
for
(
int
k
=
0
;
k
<
8
;
++
k
)
printf
(
"%g %g %g
\n
"
,
xyz
[
k
][
0
],
xyz
[
k
][
1
],
xyz
[
k
][
2
]);
}
int
main
(
int
argc
,
char
*
argv
[])
void
CubeTransform
::
LoadMapping
(
std
::
string
mapfile
,
int
dims
)
{
if
(
argc
<
3
)
{
printf
(
"Usage: %s mapping_file mesh_file [output=out.pmf] [dims=3]
\n
"
,
argv
[
0
]);
printf
(
"mapping file contains 8 (for 3D) 4 (for 2D) 2 (for 1D) coordinates mapping each corner of the cube:
\n
"
);
for
(
int
k
=
0
;
k
<
8
;
++
k
)
printf
(
"x%d y%d z%d
\n
"
,
k
,
k
,
k
);
printf
(
" (6)----------(7)
\n
"
);
printf
(
" /| /|
\n
"
);
printf
(
" / | / |
\n
"
);
printf
(
" / | / |
\n
"
);
printf
(
" (4)----------(5) |
\n
"
);
printf
(
" | | | |
\n
"
);
printf
(
" | | | |
\n
"
);
printf
(
" | (2)-------|--(3)
\n
"
);
printf
(
" | / | /
\n
"
);
printf
(
" | / | /
\n
"
);
printf
(
" |/ |/
\n
"
);
printf
(
"z (0)----------(1)
\n
"
);
printf
(
"^ y
\n
"
);
printf
(
"|/
\n
"
);
printf
(
"*-->x
\n
"
);
return
-
1
;
}
std
::
string
mapfile
(
argv
[
1
]);
std
::
string
infile
(
argv
[
2
]);
std
::
string
outfile
=
"out.pmf"
;
int
dims
=
3
;
if
(
argc
>
3
)
outfile
=
std
::
string
(
argv
[
3
]);
if
(
argc
>
4
)
dims
=
atoi
(
argv
[
4
]);
std
::
fstream
map
(
mapfile
.
c_str
(),
std
::
ios
::
in
);
double
xyz
[
8
][
3
]
=
{{
0
,
0
,
0
},{
1
,
0
,
0
},{
0
,
1
,
0
},{
1
,
1
,
0
},{
0
,
0
,
1
},{
1
,
0
,
1
},{
0
,
1
,
1
},{
1
,
1
,
1
}};
for
(
int
k
=
0
;
k
<
(
1
<<
dims
);
++
k
)
for
(
int
d
=
0
;
d
<
dims
;
++
d
)
map
>>
xyz
[
k
][
d
];
map
.
close
();
//complete definition to other dims
for
(
int
k
=
(
1
<<
dims
);
k
<
8
;
k
+=
(
1
<<
dims
))
{
for
(
int
q
=
0
;
q
<
(
1
<<
dims
);
++
q
)
...
...
@@ -59,12 +30,27 @@ int main(int argc, char *argv[])
xyz
[
k
+
q
][
d
]
=
xyz
[
k
-
(
1
<<
dims
)
+
q
][
d
];
}
for
(
int
k
=
0
;
k
<
8
;
++
k
)
printf
(
"%g %g %g
\n
"
,
xyz
[
k
][
0
],
xyz
[
k
][
1
],
xyz
[
k
][
2
]);
Mesh
mesh
;
mesh
.
Load
(
infile
);
}
void
CubeTransform
::
SetMapping
(
double
*
map
,
int
dims
)
{
int
q
=
0
;
for
(
int
k
=
0
;
k
<
(
1
<<
dims
);
++
k
)
for
(
int
d
=
0
;
d
<
dims
;
++
d
)
xyz
[
k
][
d
]
=
map
[
q
++
];
//complete definition to other dims
for
(
int
k
=
(
1
<<
dims
);
k
<
8
;
k
+=
(
1
<<
dims
))
{
for
(
int
q
=
0
;
q
<
(
1
<<
dims
);
++
q
)
for
(
int
d
=
0
;
d
<
dims
;
++
d
)
xyz
[
k
+
q
][
d
]
=
xyz
[
k
-
(
1
<<
dims
)
+
q
][
d
];
}
}
void
CubeTransform
::
Transform
()
{
double
cmin
[
3
],
cmax
[
3
];
for
(
int
d
=
0
;
d
<
3
;
++
d
)
{
...
...
@@ -90,7 +76,7 @@ int main(int argc, char *argv[])
double
c
[
3
]
=
{
0
,
0
,
0
};
for
(
int
d
=
0
;
d
<
3
;
++
d
)
c
[
d
]
=
(
n
->
Coords
()[
d
]
-
cmin
[
d
])
/
(
cmax
[
d
]
-
cmin
[
d
]);
for
(
int
d
=
0
;
d
<
dims
;
++
d
)
for
(
int
d
=
0
;
d
<
3
;
++
d
)
n
->
Coords
()[
d
]
=
((
xyz
[
0
][
d
]
*
(
1
-
c
[
0
])
+
xyz
[
1
][
d
]
*
c
[
0
])
*
(
1
-
c
[
1
])
+
(
xyz
[
2
][
d
]
*
(
1
-
c
[
0
])
+
xyz
[
3
][
d
]
*
c
[
0
])
*
c
[
1
])
*
(
1
-
c
[
2
])
+
((
xyz
[
4
][
d
]
*
(
1
-
c
[
0
])
+
xyz
[
5
][
d
]
*
c
[
0
])
*
(
1
-
c
[
1
])
+
...
...
@@ -116,14 +102,5 @@ int main(int argc, char *argv[])
for
(
int
d
=
0
;
d
<
3
;
++
d
)
printf
(
"%d %g:%g
\n
"
,
d
,
cmin
[
d
],
cmax
[
d
]);
Storage
::
bulk_array
name
=
mesh
->
self
()
->
BulkArray
(
mesh
.
CreateTag
(
"GRIDNAME"
,
DATA_BULK
,
MESH
,
NONE
));
name
.
insert
(
name
.
end
(),
mesh_name
.
begin
(),
mesh_name
.
end
());
printf
(
"I'm ready!
\n
"
);
mesh
.
Save
(
outfile
);
return
0
;
}
}
Examples/GridTools/cube_transform.h
0 → 100644
View file @
cc18b29b
#ifndef _CUBE_TRANSFORM_H
#define _CUBE_TRANSFORM_H
#include "inmost.h"
using
namespace
INMOST
;
class
CubeTransform
{
Mesh
&
mesh
;
double
xyz
[
8
][
3
];
public:
CubeTransform
(
Mesh
&
m
)
:
mesh
(
m
),
xyz
({{
0
,
0
,
0
},{
1
,
0
,
0
},{
0
,
1
,
0
},{
1
,
1
,
0
},{
0
,
0
,
1
},{
1
,
0
,
1
},{
0
,
1
,
1
},{
1
,
1
,
1
}})
{}
CubeTransform
(
const
CubeTransform
&
b
)
:
mesh
(
b
.
mesh
)
{}
CubeTransform
&
operator
=
(
CubeTransform
const
&
b
)
{
mesh
=
b
.
mesh
;
return
*
this
;}
~
CubeTransform
()
{}
/// Set mapping.
/// @param map Array map should contain (2^dims)*dims integers.
/// @param dims Number of dimensions.
void
SetMapping
(
double
*
map
,
int
dims
=
3
);
/// Load mapping from file.
/// @param file File name.
/// @param dims Number of dimensions.
void
LoadMapping
(
std
::
string
file
,
int
dims
=
3
);
/// Output mapping to std::cout.
void
PrintMapping
();
///Perform transformation.
void
Transform
();
};
#endif //_CUBE_TRANSFORM_H
Examples/GridTools/fix_faults.cpp
View file @
cc18b29b
#include "inmost.h"
#include "fix_faults.h"
using
namespace
INMOST
;
//todo: non-flat faces
...
...
@@ -448,21 +448,9 @@ struct coords
int
main
(
int
argc
,
char
**
argv
)
void
FixFaults
::
FixMeshFaults
(
)
{
if
(
argc
<
2
)
{
std
::
cout
<<
"Usage: "
<<
argv
[
0
]
<<
" mesh [mesh_out=grid.vtk]"
<<
std
::
endl
;
return
-
1
;
}
std
::
string
grid_out
=
"grid.vtk"
;
if
(
argc
>
2
)
grid_out
=
std
::
string
(
argv
[
2
]);
Mesh
m
;
m
.
Load
(
argv
[
1
]);
std
::
cout
<<
"Start:"
<<
std
::
endl
;
std
::
cout
<<
"Cells: "
<<
m
.
NumberOfCells
()
<<
std
::
endl
;
std
::
cout
<<
"Faces: "
<<
m
.
NumberOfFaces
()
<<
std
::
endl
;
...
...
@@ -919,6 +907,9 @@ int main(int argc, char ** argv)
}
}
m
.
EndModification
();
m
.
DeleteTag
(
new_points
);
m
.
DeleteTag
(
new_edges
);
std
::
cout
<<
"time to split: "
<<
Timer
()
-
tt
<<
std
::
endl
;
...
...
@@ -927,7 +918,4 @@ int main(int argc, char ** argv)
std
::
cout
<<
"Cells: "
<<
m
.
NumberOfCells
()
<<
std
::
endl
;
std
::
cout
<<
"Faces: "
<<
m
.
NumberOfFaces
()
<<
std
::
endl
;
std
::
cout
<<
"Edges: "
<<
m
.
NumberOfEdges
()
<<
std
::
endl
;
m
.
Save
(
grid_out
);
return
0
;
}
\ No newline at end of file
}
Examples/GridTools/fix_faults.h
0 → 100644
View file @
cc18b29b
#ifndef _FIX_FAULTS_H
#define _FIX_FAULTS_H
#include "inmost.h"
using
namespace
INMOST
;
class
FixFaults
{
Mesh
&
m
;
public:
FixFaults
(
Mesh
&
m
)
:
m
(
m
)
{}
FixFaults
(
const
FixFaults
&
b
)
:
m
(
b
.
m
)
{}
FixFaults
&
operator
=
(
FixFaults
const
&
b
)
{
m
=
b
.
m
;
return
*
this
;}
~
FixFaults
()
{}
void
FixMeshFaults
();
};
#endif //_FIX_FAULTS_H
Examples/GridTools/fracture.cpp
View file @
cc18b29b
...
...
@@ -72,7 +72,7 @@ void Fracture::FaceCenter(Face f, INMOST_DATA_REAL_TYPE cnt[3]) const
else
f
->
Centroid
(
cnt
);
}
void
Fracture
::
Open
(
Tag
aperture
,
bool
fill_fracture
)
void
Fracture
::
Open
(
Tag
aperture
,
bool
fill_fracture
,
double
gap_multiplier
)
{
fracture_aperture
=
aperture
;
m
->
BeginModification
();
...
...
@@ -82,7 +82,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture)
m
->
self
()
->
Integer
(
m
->
CreateTag
(
"FRACTURE_MARKER"
,
DATA_INTEGER
,
MESH
,
NONE
,
1
))
=
fracture_marker
;
//break up all the faces with fractures into volumes?
fracture_volume
=
m
->
CreateTag
(
"VOLUME_FRACTURE"
,
DATA_REAL
,
CELL
,
CELL
,
1
);
const
Storage
::
real
mult
=
0.95
;
const
Storage
::
real
mult
=
gap_multiplier
;
Tag
indexes
=
m
->
CreateTag
(
"FRAC_TEMP_INDEXES"
,
DATA_INTEGER
,
CELL
,
NONE
,
2
);
//associate a node or a pair of nodes to each non-fracture edge that ends with a fracture node
Tag
cell2node
=
m
->
CreateTag
(
"CELL2NODE"
,
DATA_REFERENCE
,
CELL
,
NONE
);
...
...
@@ -333,7 +333,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture)
}
}
}
/*
std
::
cout
<<
"Copying following real-value tags from fracture faces to fracture cells:"
<<
std
::
endl
;
for
(
std
::
vector
<
Tag
>::
iterator
q
=
transfer_real_tags
.
begin
();
q
!=
transfer_real_tags
.
end
();
++
q
)
std
::
cout
<<
q
->
GetTagName
()
<<
std
::
endl
;
...
...
@@ -341,7 +341,7 @@ void Fracture::Open(Tag aperture, bool fill_fracture)
std
::
cout
<<
"Copying following integer-value tags from fracture faces to fracture cells:"
<<
std
::
endl
;
for
(
std
::
vector
<
Tag
>::
iterator
q
=
transfer_integer_tags
.
begin
();
q
!=
transfer_integer_tags
.
end
();
++
q
)
std
::
cout
<<
q
->
GetTagName
()
<<
std
::
endl
;
*/
//now create fracture-fracture control volumes, add fracture-matrix faces to matrix cells
...
...
Examples/GridTools/fracture.h
View file @
cc18b29b
...
...
@@ -15,8 +15,14 @@ class Fracture
MarkerType
multiple_fracture_joints
;
MarkerType
matrix_fracture
;
public:
/// Marker indicates that the cell is the fracture cell.
/// For face indicates that the face appears between two
/// fracture cells.
MarkerType
isFracture
()
const
{
return
fracture_marker
;}
/// Marker indicates that there additional adjacent fracture cells
/// and multiple fracture cells meet at one of the edges of the face.
MarkerType
MultipleJoints
()
const
{
return
multiple_fracture_joints
;}
/// Marker indicates that face appears between matrix and fracture.
MarkerType
MatrixFracture
()
const
{
return
matrix_fracture
;}
///Default constructor, defines empty set.
Fracture
(
Mesh
&
_m
)
:
m
(
&
_m
)
{}
...
...
@@ -24,12 +30,23 @@ public:
Fracture
(
const
Fracture
&
b
)
:
m
(
b
.
m
)
{}
///Assignment operator.
Fracture
&
operator
=
(
Fracture
const
&
b
)
{
m
=
b
.
m
;
return
*
this
;}
/// Retrive cell volume that accounts for fracture aperture.
INMOST_DATA_REAL_TYPE
Volume
(
Cell
c
)
const
;
/// Retrive face area that accounts for fracture aperture.
INMOST_DATA_REAL_TYPE
Area
(
Face
f
)
const
;
/// Retrive face center that accounts for fracture aperture.
void
FaceCenter
(
Face
f
,
INMOST_DATA_REAL_TYPE
cnt
[
3
])
const
;
void
Open
(
Tag
aperture
,
bool
fill_fracture
);
/// @param aperture Defines the opening of the fracture. Geometry does not
/// necesserely follows this parameter, however the functions for Volume and
/// Area are affected by this parameter.
/// @param fill_fracture Tells whether it is needed to connect the mesh with cells
/// in the fracture region. If parameter is false the mesh is disconnected.
/// @param gap_multiplier This parameter averages the placement of the new mesh node
/// position between the original node of the fracture face and averaged centers of cells
/// that appear on the same side of fracture. If the parameter is 1, then the new node
/// is positioned at the same location as original node.
void
Open
(
Tag
aperture
,
bool
fill_fracture
,
double
gap_multiplier
=
0.95
);
};
#endif //__FRACTURE_H
Examples/GridTools/set_layers.cpp
0 → 100644
View file @
cc18b29b
#include "set_layers.h"
#define ind(r,c) ((r)*N + (c))
void
rand1d
(
double
*
arr
,
int
N
)
{
if
(
N
==
1
)
return
;
const
double
t
=
0.15
;
double
a
=
0.5
+
(
2
*
(
rand
()
*
1.0
/
RAND_MAX
)
-
1
)
*
t
;
// in 0.5+[-1,1]*t
arr
[
N
/
2
]
=
arr
[
0
]
*
a
+
arr
[
N
]
*
(
1
-
a
);
rand1d
(
arr
,
N
/
2
);
rand1d
(
arr
+
N
/
2
,
N
-
N
/
2
);
}
std
::
pair
<
int
,
double
>
layer1d
(
double
*
arr
,
int
N
,
double
z
)
{
for
(
int
k
=
0
;
k
<
N
;
++
k
)
{
if
(
z
>=
arr
[
k
]
&&
z
<=
arr
[
k
+
1
]
)
return
std
::
make_pair
(
k
,
(
z
-
arr
[
k
])
/
(
arr
[
k
+
1
]
-
arr
[
k
]));
}
return
std
::
make_pair
(
-
1
,
0.0
);
}
const
double
noind
=
std
::
numeric_limits
<
double
>::
quiet_NaN
();
void
rand2d
(
double
*
arr
,
int
N
,
int
Nl
,
int
Nr
,
int
Nb
,
int
Nt
,
double
t
)
{
//std::cout << "Nl:Nr " << Nl <<":" << Nr << " Nb:Nt " << Nb << ":" << Nt << std::endl;
if
(
Nr
-
Nl
<
2
&&
Nt
-
Nb
<
2
)
{
//std::cout << "exit" << std::endl;
return
;
}
//const double t = 0.15;
int
Nk
=
(
Nb
+
Nt
)
/
2
;
int
Nm
=
(
Nl
+
Nr
)
/
2
;
//std::cout << "Nk " << Nk << " Nm " << Nm << std::endl;
double
lb
=
arr
[
ind
(
Nb
,
Nl
)];
double
rb
=
arr
[
ind
(
Nb
,
Nr
)];
double
lt
=
arr
[
ind
(
Nt
,
Nl
)];
double
rt
=
arr
[
ind
(
Nt
,
Nr
)];
if
(
lb
!=
lb
||
rb
!=
rb
||
lt
!=
lt
||
rt
!=
rt
)
throw
-
1
;
if
(
arr
[
ind
(
Nk
,
Nl
)]
!=
arr
[
ind
(
Nk
,
Nl
)]
)
arr
[
ind
(
Nk
,
Nl
)]
=
0.5
*
(
lb
+
lt
)
+
(
2
*
(
rand
()
*
1.0
/
RAND_MAX
)
-
1
)
*
t
;
if
(
arr
[
ind
(
Nk
,
Nr
)]
!=
arr
[
ind
(
Nk
,
Nr
)]
)
arr
[
ind
(
Nk
,
Nr
)]
=
0.5
*
(
rb
+
rt
)
+
(
2
*
(
rand
()
*
1.0
/
RAND_MAX
)
-
1
)
*
t
;
if
(
arr
[
ind
(
Nb
,
Nm
)]
!=
arr
[
ind
(
Nb
,
Nm
)]
)
arr
[
ind
(
Nb
,
Nm
)]
=
0.5
*
(
lb
+
rb
)
+
(
2
*
(
rand
()
*
1.0
/
RAND_MAX
)
-
1
)
*
t
;
if
(
arr
[
ind
(
Nt
,
Nm
)]
!=
arr
[
ind
(
Nt
,
Nm
)]
)
arr
[
ind
(
Nt
,
Nm
)]
=
0.5
*
(
lt
+
rt
)
+
(
2
*
(
rand
()
*
1.0
/
RAND_MAX
)
-
1
)
*
t
;
arr
[
ind
(
Nk
,
Nm
)]
=
0.25
*
(
lb
+
rb
+
lt
+
rt
)
+
(
2
*
(
rand
()
*
1.0
/
RAND_MAX
)
-
1
)
*
t
;
rand2d
(
arr
,
N
,
Nl
,
Nm
,
Nb
,
Nk
,
t
*
0.5
);
rand2d
(
arr
,
N
,
Nm
,
Nr
,
Nb
,
Nk
,
t
*
0.5
);
rand2d
(
arr
,
N
,
Nl
,
Nm
,
Nk
,
Nt
,
t
*
0.5
);
rand2d
(
arr
,
N
,
Nm
,
Nr
,
Nk
,
Nt
,
t
*
0.5
);
}
void
init2d
(
double
*
arr
,
int
N
,
double
mint
,
double
maxt
)
{
for
(
int
k
=
0
;
k
<
N
*
N
;
++
k
)
arr
[
k
]
=
noind
;
arr
[
ind
(
0
,
0
)]
=
mint
+
(
rand
()
*
1.0
/
RAND_MAX
)
*
(
maxt
-
mint
);
arr
[
ind
(
0
,
N
-
1
)]
=
mint
+
(
rand
()
*
1.0
/
RAND_MAX
)
*
(
maxt
-
mint
);
arr
[
ind
(
N
-
1
,
0
)]
=
mint
+
(
rand
()
*
1.0
/
RAND_MAX
)
*
(
maxt
-
mint
);
arr
[
ind
(
N
-
1
,
N
-
1
)]
=
mint
+
(
rand
()
*
1.0
/
RAND_MAX
)
*
(
maxt
-
mint
);
}
double
intrp2d
(
double
*
arr
,
int
N
,
double
x
,
double
y
)
{
int
n
=
ceil
(
x
*
(
N
-
1
));
int
m
=
ceil
(
y
*
(
N
-
1
));
if
(
n
==
0
)
n
=
1
;
if
(
m
==
0
)
m
=
1
;
double
dh
=
1.0
/
(
double
)(
N
-
1
);
double
kx
=
(
x
-
(
n
-
1
)
*
dh
)
/
dh
;
double
ky
=
(
y
-
(
m
-
1
)
*
dh
)
/
dh
;
//if( kx < 0 || kx > 1 ) std::cout << "bad kx: " << kx << " x is " << x << " n is " << n << " dh is " << dh << " N is " << N << std::endl;
//if( ky < 0 || ky > 1 ) std::cout << "bad ky: " << ky << " y is " << y << " m is " << m << " dh is " << dh << " N is " << N << std::endl;
double
lb
=
arr
[
ind
(
m
-
1
,
n
-
1
)];
double
rb
=
arr
[
ind
(
m
-
1
,
n
+
0
)];
double
lt
=
arr
[
ind
(
m
+
0
,
n
-
1
)];
double
rt
=
arr
[
ind
(
m
+
0
,
n
+
0
)];
if
(
lb
!=
lb
||
rb
!=
rb
||
lt
!=
lt
||
rt
!=
rt
)
throw
-
1
;
return
(
1
-
ky
)
*
(
lb
*
(
1
-
kx
)
+
rb
*
kx
)
+
ky
*
(
lt
*
(
1
-
kx
)
+
rt
*
kx
);
}
void
SetLayers
::
RandomLayersZ
(
int
num_layers
)
{
layers_z
.
resize
(
num_layers
+
1
);
layers_z
[
0
]
=
0
;
layers_z
[
num_layers
]
=
1
;
rand1d
(
&
layers_z
[
0
],
num_layers
);
}
void
SetLayers
::
DeformHeight
(
double
height
)
{
const
int
N
=
128
;
std
::
vector
<
double
>
map
(
N
*
N
,
0.0
);
init2d
(
&
map
[
0
],
N
,
0.0
,
height
);
rand2d
(
&
map
[
0
],
N
,
0
,
N
-
1
,
0
,
N
-
1
,
height
*
0.5
);
for
(
int
k
=
0
;
k
<
N
*
N
;
++
k
)
if
(
map
[
k
]
!=
map
[
k
]
)
std
::
cout
<<
"NaN at row "
<<
k
/
N
<<
" col "
<<
k
%
N
<<
std
::
endl
;
double
cmin
[
3
],
cmax
[
3
];
for
(
int
d
=
0
;
d
<
3
;
++
d
)
{
cmin
[
d
]
=
1.0e+20
;
cmax
[
d
]
=-
1.0e+20
;
}
for
(
Mesh
::
iteratorNode
n
=
mesh
.
BeginNode
();
n
!=
mesh
.
EndNode
();
++
n
)
{
for
(
int
d
=
0
;
d
<
3
;
++
d
)
{
if
(
cmin
[
d
]
>
n
->
Coords
()[
d
]
)
cmin
[
d
]
=
n
->
Coords
()[
d
];
if
(
cmax
[
d
]
<
n
->
Coords
()[
d
]
)
cmax
[
d
]
=
n
->
Coords
()[
d
];
}
}
for
(
Mesh
::
iteratorNode
n
=
mesh
.
BeginNode
();
n
!=
mesh
.
EndNode
();
++
n
)
{
double
c
[
3
]
=
{
0
,
0
,
0
};
for
(
int
d
=
0
;
d
<
3
;
++
d
)
c
[
d
]
=
(
n
->
Coords
()[
d
]
-
cmin
[
d
])
/
(
cmax
[
d
]
-
cmin
[
d
]);
n
->
Coords
()[
2
]
+=
intrp2d
(
&
map
[
0
],
N
,
c
[
0
],
c
[
1
]);
}
}
void
SetLayers
::
DeformLayers
(
double
coef
)
{
const
int
N
=
128
;
std
::
vector
<
std
::
vector
<
double
>
>
map_layer
(
layers_z
.
size
());
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for
(
int
k
=
0
;
k
<
map_layer
.
size
();
++
k
)
{
map_layer
[
k
].
resize
(
N
*
N
);
init2d
(
&
map_layer
[
k
][
0
],
N
,
1
-
coef
,
1
);
rand2d
(
&
map_layer
[
k
][
0
],
N
,
0
,
N
-
1
,
0
,
N
-
1
,
coef
);
for
(
int
kk
=
0
;
kk
<
N
*
N
;
++
kk
)
if
(
map_layer
[
k
][
kk
]
<
0.0
)
map_layer
[
k
][
kk
]
=
0
;
}
double
cmin
[
3
],
cmax
[
3
];
for
(
int
d
=
0
;
d
<
3
;
++
d
)
{
cmin
[
d
]
=
1.0e+20
;
cmax
[
d
]
=-
1.0e+20
;
}
for
(
Mesh
::
iteratorNode
n
=
mesh
.
BeginNode
();
n
!=
mesh
.
EndNode
();
++
n
)
{
for
(
int
d
=
0
;
d
<
3
;
++
d
)
{
if
(
cmin
[
d
]
>
n
->
Coords
()[
d
]
)
cmin
[
d
]
=
n
->
Coords
()[
d
];
if
(
cmax
[
d
]
<
n
->
Coords
()[
d
]
)
cmax
[
d
]
=
n
->
Coords
()[
d
];
}
}
TagInteger
layer_tag
=
mesh
.
CreateTag
(
"LAYER"
,
DATA_INTEGER
,
CELL
,
NONE
,
1
);
TagReal
coef_tag
=
mesh
.
CreateTag
(
"LAYER_COEF"
,
DATA_REAL
,
CELL
,
NONE
,
1
);