Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
fluidparticles
MigFlow
Commits
eb12c274
Commit
eb12c274
authored
Nov 14, 2018
by
Jonathan Lambrechts
Browse files
Merge branch 'oneFluid' into weakBnd
parents
1ad4ac09
f7b61685
Pipeline
#4584
passed with stage
in 1 minute and 48 seconds
Changes
14
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
fluid/CMakeLists.txt
View file @
eb12c274
...
...
@@ -27,7 +27,6 @@ set(SRC
hxt_linear_system.c
hxt_linear_system_lu.c
hxt_message.c
fluid_problem_io.c
)
if
(
PETSC_FOUND
)
...
...
fluid/fluid_problem.c
View file @
eb12c274
...
...
@@ -1141,6 +1141,7 @@ void fluid_problem_free(FluidProblem *problem)
free
(
problem
->
solution
);
free
(
problem
->
porosity
);
free
(
problem
->
oldporosity
);
free
(
problem
->
node_volume
);
hxtLinearSystemDelete
(
&
problem
->
linear_system
);
free
(
problem
->
particle_uvw
);
free
(
problem
->
particle_element_id
);
...
...
@@ -1161,6 +1162,7 @@ void fluid_problem_free(FluidProblem *problem)
free
(
problem
->
strong_boundaries
);
mesh_tree_free
(
problem
->
mesh_tree
);
mesh_free
(
problem
->
mesh
);
free
(
problem
);
}
void
fluid_problem_adapt_gen_mesh
(
FluidProblem
*
problem
,
double
lcmax
,
double
lcmin
,
double
n_el
,
double
cmax
,
double
cmin
)
...
...
@@ -1594,3 +1596,81 @@ void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const cha
p
->
weak_boundaries
[
p
->
n_weak_boundaries
-
1
].
cb
=
f
;
p
->
weak_boundaries
[
p
->
n_weak_boundaries
-
1
].
apply
=
apply
;
}
int
fluid_problem_n_elements
(
const
FluidProblem
*
p
)
{
return
p
->
mesh
->
n_elements
;
}
int
*
fluid_problem_elements
(
const
FluidProblem
*
p
)
{
return
p
->
mesh
->
elements
;
}
double
*
fluid_problem_porosity
(
const
FluidProblem
*
p
)
{
return
p
->
porosity
;
}
double
*
fluid_problem_old_porosity
(
const
FluidProblem
*
p
)
{
return
p
->
oldporosity
;
}
int
fluid_problem_private_n_physical
(
const
FluidProblem
*
p
){
return
p
->
mesh
->
n_phys
;
}
void
fluid_problem_private_physical
(
const
FluidProblem
*
p
,
int
i
,
char
**
phys_name
,
int
*
phys_dim
,
int
*
phys_id
,
int
*
phys_n_nodes
,
int
**
phys_nodes
,
int
*
phys_n_boundaries
,
int
**
phys_boundaries
)
{
*
phys_name
=
p
->
mesh
->
phys_name
[
i
];
*
phys_n_nodes
=
p
->
mesh
->
phys_n_nodes
[
i
];
*
phys_nodes
=
p
->
mesh
->
phys_nodes
[
i
];
*
phys_dim
=
p
->
mesh
->
phys_dimension
[
i
];
*
phys_id
=
p
->
mesh
->
phys_id
[
i
];
*
phys_n_boundaries
=
p
->
mesh
->
phys_n_boundaries
[
i
];
*
phys_boundaries
=
p
->
mesh
->
phys_boundaries
[
i
];
}
void
fluid_problem_private_set_elements
(
FluidProblem
*
p
,
int
n_nodes
,
double
*
x
,
int
n_elements
,
int
*
elements
)
{
Mesh
*
m
=
malloc
(
sizeof
(
Mesh
));
m
->
n_phys
=
0
;
m
->
phys_name
=
NULL
;
m
->
phys_id
=
NULL
;
m
->
phys_dimension
=
NULL
;
m
->
phys_n_nodes
=
NULL
;
m
->
phys_nodes
=
NULL
;
m
->
phys_n_boundaries
=
NULL
;
m
->
phys_boundaries
=
NULL
;
m
->
n_nodes
=
n_nodes
;
m
->
n_elements
=
n_elements
;
m
->
x
=
malloc
(
n_nodes
*
3
*
sizeof
(
double
));
memcpy
(
m
->
x
,
x
,
sizeof
(
double
)
*
3
*
n_nodes
);
m
->
elements
=
malloc
(
sizeof
(
int
)
*
(
DIMENSION
+
1
)
*
n_elements
);
memcpy
(
m
->
elements
,
elements
,(
DIMENSION
+
1
)
*
sizeof
(
int
)
*
n_elements
);
fluid_problem_set_mesh
(
p
,
m
);
}
void
fluid_problem_private_add_physical
(
FluidProblem
*
p
,
int
tag
,
int
dim
,
const
char
*
name
,
int
n_nodes
,
int
*
nodes
,
int
n_bnd
,
int
*
bnd
)
{
Mesh
*
m
=
p
->
mesh
;
int
i
=
m
->
n_phys
;
m
->
n_phys
+=
1
;
m
->
phys_name
=
realloc
(
m
->
phys_name
,
m
->
n_phys
*
sizeof
(
char
*
));
m
->
phys_name
[
i
]
=
strdup
(
name
);
m
->
phys_id
=
realloc
(
m
->
phys_id
,
m
->
n_phys
*
sizeof
(
int
));
m
->
phys_id
[
i
]
=
tag
;
m
->
phys_dimension
=
realloc
(
m
->
phys_dimension
,
m
->
n_phys
*
sizeof
(
int
));
m
->
phys_dimension
[
i
]
=
dim
;
m
->
phys_n_nodes
=
realloc
(
m
->
phys_n_nodes
,
m
->
n_phys
*
sizeof
(
int
));
m
->
phys_n_nodes
[
i
]
=
n_nodes
;
m
->
phys_nodes
=
realloc
(
m
->
phys_nodes
,
m
->
n_phys
*
sizeof
(
int
*
));
m
->
phys_nodes
[
i
]
=
malloc
(
n_nodes
*
sizeof
(
int
));
memcpy
(
m
->
phys_nodes
[
i
],
nodes
,
n_nodes
*
sizeof
(
int
));
m
->
phys_n_boundaries
=
realloc
(
m
->
phys_n_boundaries
,
m
->
n_phys
*
sizeof
(
int
));
m
->
phys_n_boundaries
[
i
]
=
n_bnd
;
m
->
phys_boundaries
=
realloc
(
m
->
phys_boundaries
,
m
->
n_phys
*
sizeof
(
int
*
));
m
->
phys_boundaries
[
i
]
=
malloc
(
n_bnd
*
sizeof
(
int
)
*
2
);
memcpy
(
m
->
phys_boundaries
[
i
],
bnd
,
n_bnd
*
sizeof
(
int
)
*
2
);
}
fluid/fluid_problem.h
View file @
eb12c274
...
...
@@ -85,9 +85,6 @@ struct FluidProblem {
int
*
particle_element_id
;
};
int
fluid_problem_export
(
const
FluidProblem
*
problem
,
const
char
*
output_dir
,
double
t
,
int
iter
);
int
fluid_problem_export_vtk
(
const
FluidProblem
*
problem
,
const
char
*
output_dir
,
double
t
,
int
iter
);
int
fluid_problem_import_vtk
(
FluidProblem
*
problem
,
const
char
*
filename
);
// complete force on the particle (including gravity)
void
fluid_problem_compute_node_particle_force
(
FluidProblem
*
problem
,
double
dt
,
double
*
particle_force
);
int
fluid_problem_implicit_euler
(
FluidProblem
*
problem
,
double
dt
,
double
newton_atol
,
double
newton_rtol
,
int
newton_max_it
,
int
reduced_gravity
,
double
stab_param
);
...
...
@@ -105,4 +102,13 @@ double *fluid_problem_coordinates(const FluidProblem *p);
int
fluid_problem_n_nodes
(
const
FluidProblem
*
p
);
void
fluid_problem_set_weak_boundary
(
FluidProblem
*
p
,
const
char
*
tag
,
const
char
*
bndtype
,
StrongBoundaryCallback
*
apply
);
void
fluid_problem_set_strong_boundary
(
FluidProblem
*
p
,
const
char
*
tag
,
int
field
,
int
row
,
StrongBoundaryCallback
*
apply
);
int
fluid_problem_n_elements
(
const
FluidProblem
*
p
);
int
*
fluid_problem_elements
(
const
FluidProblem
*
p
);
double
*
fluid_problem_old_porosity
(
const
FluidProblem
*
p
);
double
*
fluid_problem_porosity
(
const
FluidProblem
*
p
);
int
fluid_problem_private_n_physical
(
const
FluidProblem
*
p
);
void
fluid_problem_private_physical
(
const
FluidProblem
*
p
,
int
i
,
char
**
phys_name
,
int
*
phys_dim
,
int
*
phys_id
,
int
*
phys_n_nodes
,
int
**
phys_nodes
,
int
*
phys_n_boundaries
,
int
**
phys_boundaries
);
void
fluid_problem_private_set_elements
(
FluidProblem
*
p
,
int
n_nodes
,
double
*
x
,
int
n_elements
,
int
*
elements
);
void
fluid_problem_private_add_physical
(
FluidProblem
*
p
,
int
tag
,
int
dim
,
const
char
*
name
,
int
n_nodes
,
int
*
nodes
,
int
n_bnd
,
int
*
bnd
);
#endif
fluid/fluid_problem_io.c
deleted
100644 → 0
View file @
1ad4ac09
/*
* MigFlow - Copyright (C) <2010-2018>
* <Universite catholique de Louvain (UCL), Belgium
* Universite de Montpellier, France>
*
* List of the contributors to the development of MigFlow: see AUTHORS file.
* Description and complete License: see LICENSE file.
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation, either version
* 3 of the License, or (at your option) any later version.
*
* This program 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 Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program (see COPYING and COPYING.LESSER files). If not,
* see <http://www.gnu.org/licenses/>.
*/
#include "fluid_problem.h"
#include "mbxml.h"
#include <string.h>
int
fluid_problem_import_vtk
(
FluidProblem
*
problem
,
const
char
*
filename
)
{
const
int
n_fields
=
fluid_problem_n_fields
(
problem
);
MbXmlReader
*
r
=
mb_xml_reader_create
(
filename
);
MbXmlElement
fileelement
,
gridelement
,
pieceelement
;
mb_xml_read_element
(
r
,
&
fileelement
);
mb_xml_read_element
(
r
,
&
gridelement
);
mb_xml_read_element
(
r
,
&
pieceelement
);
if
(
problem
->
mesh
)
mesh_free
(
problem
->
mesh
);
problem
->
mesh
=
malloc
(
sizeof
(
Mesh
));
Mesh
*
m
=
problem
->
mesh
;
m
->
n_nodes
=
-
1
,
m
->
n_elements
=
-
1
;
for
(
size_t
i
=
0
;
i
<
pieceelement
.
n_attr
;
++
i
)
{
if
(
!
strcasecmp
(
"NumberOfPoints"
,
pieceelement
.
attr
[
i
].
name
))
sscanf
(
pieceelement
.
attr
[
i
].
value
,
"%i"
,
&
m
->
n_nodes
);
if
(
!
strcasecmp
(
"NumberOfCells"
,
pieceelement
.
attr
[
i
].
name
))
sscanf
(
pieceelement
.
attr
[
i
].
value
,
"%i"
,
&
m
->
n_elements
);
}
if
(
m
->
n_nodes
<
0
||
m
->
n_elements
<
0
){
printf
(
"NumberOfPoints or NumberOfCells not specified in vtk piece"
);
}
MBVtkDataArray
a
;
MbXmlElement
e
;
while
(
mb_xml_read_element
(
r
,
&
e
)
!=
1
)
{
if
(
!
strcasecmp
(
e
.
name
,
"Points"
)){
mb_read_vtk_data_array
(
r
,
&
a
);
m
->
x
=
mb_vtk_data_array_to_double
(
&
a
,
m
->
n_nodes
,
3
);
mb_vtk_data_array_destroy
(
&
a
);
}
else
if
(
!
strcasecmp
(
e
.
name
,
"Cells"
))
{
int
ncon
=
-
1
;
while
(
mb_read_vtk_data_array
(
r
,
&
a
)
!=
0
){
if
(
!
strcasecmp
(
a
.
name
,
"Connectivity"
))
{
m
->
elements
=
mb_vtk_data_array_to_int
(
&
a
,
m
->
n_elements
,
DIMENSION
+
1
);
mb_vtk_data_array_destroy
(
&
a
);
}
else
if
(
!
strcasecmp
(
a
.
name
,
"Offsets"
)){
int
*
offsets
=
mb_vtk_data_array_to_int
(
&
a
,
m
->
n_elements
,
1
);
ncon
=
offsets
[
m
->
n_elements
-
1
];
free
(
offsets
);
mb_vtk_data_array_destroy
(
&
a
);
}
else
if
(
!
strcasecmp
(
a
.
name
,
"Types"
)){
mb_vtk_data_array_destroy
(
&
a
);
}
else
{
printf
(
"unknown data array
\"
%s
\"
in vtk piece
\n
"
,
a
.
name
);
exit
(
1
);
}
}
}
else
if
(
!
strcasecmp
(
e
.
name
,
"PointData"
))
{
if
(
problem
->
solution
){
free
(
problem
->
solution
);
}
problem
->
solution
=
(
double
*
)
malloc
(
sizeof
(
double
)
*
n_fields
*
m
->
n_nodes
);
while
(
mb_read_vtk_data_array
(
r
,
&
a
)
!=
0
){
if
(
!
strcasecmp
(
a
.
name
,
"Porosity"
))
{
free
(
problem
->
porosity
);
problem
->
porosity
=
mb_vtk_data_array_to_double
(
&
a
,
m
->
n_nodes
,
1
);
}
else
if
(
!
strcasecmp
(
a
.
name
,
"OldPorosity"
))
{
free
(
problem
->
oldporosity
);
problem
->
oldporosity
=
mb_vtk_data_array_to_double
(
&
a
,
m
->
n_nodes
,
1
);
}
else
if
(
!
strcasecmp
(
a
.
name
,
"Pressure"
))
{
double
*
p
=
mb_vtk_data_array_to_double
(
&
a
,
m
->
n_nodes
,
1
);
for
(
int
i
=
0
;
i
<
m
->
n_nodes
;
++
i
)
problem
->
solution
[
n_fields
*
i
+
DIMENSION
]
=
p
[
i
];
free
(
p
);
}
else
if
(
!
strcasecmp
(
a
.
name
,
"Velocity"
))
{
double
*
v
=
mb_vtk_data_array_to_double
(
&
a
,
m
->
n_nodes
,
3
);
for
(
int
i
=
0
;
i
<
m
->
n_nodes
;
++
i
)
for
(
int
d
=
0
;
d
<
DIMENSION
;
++
d
)
problem
->
solution
[
n_fields
*
i
+
d
]
=
v
[
i
*
3
+
d
];
free
(
v
);
}
else
if
(
!
strcasecmp
(
a
.
name
,
"Concentration"
))
{
double
*
c
=
mb_vtk_data_array_to_double
(
&
a
,
m
->
n_nodes
,
1
);
for
(
int
i
=
0
;
i
<
m
->
n_nodes
;
++
i
)
problem
->
solution
[
n_fields
*
i
+
DIMENSION
+
1
]
=
c
[
i
];
free
(
c
);
}
mb_vtk_data_array_destroy
(
&
a
);
}
}
else
{
printf
(
"unknown type in piece :
\"
%s
\"\n
"
,
e
.
name
);
exit
(
1
);
}
mb_xml_end_element
(
r
,
&
e
);
}
mb_xml_end_element
(
r
,
&
pieceelement
);
MbXmlElement
fielddataelement
;
mb_xml_read_element
(
r
,
&
fielddataelement
);
int
n_phys
=
0
;
int
n_phys_alloc
=
10
;
MBVtkDataArray
*
phys
=
(
MBVtkDataArray
*
)
malloc
(
sizeof
(
MBVtkDataArray
)
*
10
);
int
n_boundaries
=
0
;
int
n_boundaries_alloc
=
10
;
MBVtkDataArray
*
boundaries
=
(
MBVtkDataArray
*
)
malloc
(
sizeof
(
MBVtkDataArray
)
*
10
);
while
(
mb_read_vtk_data_array
(
r
,
&
a
)
!=
0
){
if
(
!
strncasecmp
(
a
.
name
,
"Physical "
,
9
))
{
if
(
n_phys
==
n_phys_alloc
)
{
n_phys_alloc
*=
2
;
phys
=
(
MBVtkDataArray
*
)
realloc
(
phys
,
n_phys_alloc
*
sizeof
(
MBVtkDataArray
));
}
phys
[
n_phys
++
]
=
a
;
}
else
if
(
!
strncasecmp
(
a
.
name
,
"Boundary "
,
9
))
{
if
(
n_boundaries
==
n_boundaries_alloc
)
{
n_boundaries_alloc
*=
2
;
boundaries
=
(
MBVtkDataArray
*
)
realloc
(
boundaries
,
n_boundaries_alloc
*
sizeof
(
MBVtkDataArray
));
}
boundaries
[
n_boundaries
++
]
=
a
;
}
else
{
printf
(
"unknown data array
\"
%s
\"
in vtk piece
\n
"
,
a
.
name
);
exit
(
1
);
}
}
m
->
phys_id
=
malloc
(
n_phys
*
sizeof
(
int
));
m
->
phys_dimension
=
(
int
*
)
malloc
(
n_phys
*
sizeof
(
int
));
m
->
phys_n_nodes
=
(
int
*
)
malloc
(
n_phys
*
sizeof
(
int
));
m
->
phys_name
=
(
char
**
)
malloc
(
n_phys
*
sizeof
(
char
*
));
m
->
phys_nodes
=
(
int
**
)
malloc
(
n_phys
*
sizeof
(
int
*
));
m
->
phys_n_boundaries
=
(
int
*
)
malloc
(
n_phys
*
sizeof
(
int
));
m
->
phys_boundaries
=
(
int
**
)
malloc
(
n_phys
*
sizeof
(
int
*
));
if
(
n_phys
!=
n_boundaries
)
{
printf
(
"import vtk errors : number of physical and boundaries data does not match"
);
}
m
->
n_phys
=
n_phys
;
for
(
int
i
=
0
;
i
<
n_phys
;
++
i
)
{
MBVtkDataArray
a
=
phys
[
i
];
int
p
,
n
;
sscanf
(
a
.
name
+
9
,
"%i %i%n"
,
&
m
->
phys_dimension
[
i
],
&
m
->
phys_id
[
i
],
&
p
);
m
->
phys_nodes
[
i
]
=
mb_vtk_data_array_to_int
(
&
phys
[
i
],
-
1
,
-
1
);
m
->
phys_n_nodes
[
i
]
=
phys
[
i
].
number_of_tuples
;
m
->
phys_name
[
i
]
=
strdup
(
a
.
name
+
p
+
10
);
mb_vtk_data_array_destroy
(
&
phys
[
i
]);
}
free
(
phys
);
for
(
int
i
=
0
;
i
<
n_boundaries
;
++
i
)
{
MBVtkDataArray
a
=
boundaries
[
i
];
int
p
,
n
;
sscanf
(
a
.
name
+
9
,
"%i %i%n"
,
&
m
->
phys_dimension
[
i
],
&
m
->
phys_id
[
i
],
&
p
);
m
->
phys_boundaries
[
i
]
=
mb_vtk_data_array_to_int
(
&
boundaries
[
i
],
-
1
,
-
1
);
m
->
phys_n_boundaries
[
i
]
=
boundaries
[
i
].
number_of_tuples
;
mb_vtk_data_array_destroy
(
&
boundaries
[
i
]);
}
free
(
boundaries
);
mb_xml_end_element
(
r
,
&
fielddataelement
);
mb_xml_end_element
(
r
,
&
gridelement
);
mb_xml_end_element
(
r
,
&
fileelement
);
mb_xml_reader_destroy
(
r
);
fluid_problem_after_import
(
problem
);
return
0
;
}
int
fluid_problem_export_vtk
(
const
FluidProblem
*
problem
,
const
char
*
output_dir
,
double
t
,
int
iter
)
{
const
int
n_fields
=
fluid_problem_n_fields
(
problem
);
Mesh
*
mesh
=
problem
->
mesh
;
char
file_name
[
1024
];
sprintf
(
file_name
,
"%s/fluid_%05i.vtu"
,
output_dir
,
iter
);
FILE
*
f
=
fopen
(
file_name
,
"w"
);
if
(
!
f
){
printf
(
"Cannot open file
\"
%s
\"
for writing.
\n
"
,
file_name
);
return
-
1
;
}
fprintf
(
f
,
"<VTKFile type=
\"
UnstructuredGrid
\"
format=
\"
0.1
\"
>
\n
"
);
fprintf
(
f
,
"<UnstructuredGrid>
\n
"
);
fprintf
(
f
,
"<Piece NumberOfPoints=
\"
%i
\"
NumberOfCells=
\"
%i
\"
>
\n
"
,
mesh
->
n_nodes
,
mesh
->
n_elements
);
fprintf
(
f
,
"<Points>
\n
"
);
fprintf
(
f
,
"<DataArray NumberOfComponents=
\"
3
\"
type=
\"
Float64
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
3
*
mesh
->
n_nodes
;
++
i
)
{
fprintf
(
f
,
"%.16e "
,
mesh
->
x
[
i
]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"</Points>
\n
"
);
fprintf
(
f
,
"<Cells>
\n
"
);
fprintf
(
f
,
"<DataArray type=
\"
Int32
\"
Name=
\"
connectivity
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_elements
*
(
DIMENSION
+
1
);
++
i
)
{
fprintf
(
f
,
"%i "
,
mesh
->
elements
[
i
]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"<DataArray type=
\"
Int32
\"
Name=
\"
offsets
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_elements
;
++
i
)
{
fprintf
(
f
,
"%i "
,
(
i
+
1
)
*
(
DIMENSION
+
1
));
}
#if DIMENSION == 2
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"<DataArray type=
\"
Int32
\"
Name=
\"
types
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_elements
;
++
i
)
{
fprintf
(
f
,
"%i "
,
5
);
// 10 for 3d
}
#else
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"<DataArray type=
\"
Int32
\"
Name=
\"
types
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_elements
;
++
i
)
{
fprintf
(
f
,
"%i "
,
10
);
// 10 for 3d
}
#endif
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"</Cells>
\n
"
);
fprintf
(
f
,
"<PointData Scalars=
\"
Porosity OldPorosity Pressure
\"
Vectors=
\"
Velocity
\"
>
\n
"
);
fprintf
(
f
,
"<DataArray Name=
\"
Porosity
\"
NumberOfComponents =
\"
1
\"
type=
\"
Float64
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_nodes
;
++
i
)
{
fprintf
(
f
,
"%.16e "
,
problem
->
porosity
[
i
]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"<DataArray Name=
\"
OldPorosity
\"
NumberOfComponents =
\"
1
\"
type=
\"
Float64
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_nodes
;
++
i
)
{
fprintf
(
f
,
"%.16e "
,
problem
->
oldporosity
[
i
]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"<DataArray Name=
\"
Pressure
\"
NumberOfComponents =
\"
1
\"
type=
\"
Float64
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_nodes
;
++
i
)
{
fprintf
(
f
,
"%.16e "
,
problem
->
solution
[
i
*
n_fields
+
(
DIMENSION
)]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
if
(
problem
->
n_fluids
==
2
){
fprintf
(
f
,
"<DataArray Name=
\"
Concentration
\"
NumberOfComponents =
\"
1
\"
type=
\"
Float64
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_nodes
;
++
i
)
{
fprintf
(
f
,
"%.16e "
,
problem
->
solution
[
i
*
n_fields
+
(
DIMENSION
+
1
)]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
}
fprintf
(
f
,
"<DataArray Name=
\"
Velocity
\"
NumberOfComponents =
\"
3
\"
type=
\"
Float64
\"
format=
\"
ascii
\"
>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_nodes
;
++
i
)
{
for
(
int
j
=
0
;
j
<
DIMENSION
;
++
j
)
fprintf
(
f
,
"%.16e "
,
problem
->
solution
[
i
*
n_fields
+
j
]);
if
(
DIMENSION
==
2
)
fprintf
(
f
,
"0 "
);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"</PointData>
\n
"
);
fprintf
(
f
,
"</Piece>
\n
"
);
fprintf
(
f
,
"<FieldData>
\n
"
);
for
(
int
i
=
0
;
i
<
mesh
->
n_phys
;
++
i
)
{
fprintf
(
f
,
"<DataArray Name=
\"
Physical %i %i %s
\"
NumberOfTuples=
\"
%i
\"
NumberOfComponents=
\"
1
\"
type=
\"
Int32
\"
format=
\"
ascii
\"
>
\n
"
,
mesh
->
phys_dimension
[
i
],
mesh
->
phys_id
[
i
],
mesh
->
phys_name
[
i
],
mesh
->
phys_n_nodes
[
i
]);
for
(
int
j
=
0
;
j
<
mesh
->
phys_n_nodes
[
i
];
++
j
)
{
fprintf
(
f
,
"%i "
,
mesh
->
phys_nodes
[
i
][
j
]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
fprintf
(
f
,
"<DataArray Name=
\"
Boundary %i %i %s
\"
NumberOfTuples=
\"
%i
\"
NumberOfComponents=
\"
2
\"
type=
\"
Int32
\"
format=
\"
ascii
\"
>
\n
"
,
mesh
->
phys_dimension
[
i
],
mesh
->
phys_id
[
i
],
mesh
->
phys_name
[
i
],
mesh
->
phys_n_boundaries
[
i
]);
for
(
int
j
=
0
;
j
<
mesh
->
phys_n_boundaries
[
i
];
++
j
)
{
fprintf
(
f
,
"%i %i "
,
mesh
->
phys_boundaries
[
i
][
j
*
2
+
0
],
mesh
->
phys_boundaries
[
i
][
j
*
2
+
1
]);
}
fprintf
(
f
,
"
\n
</DataArray>
\n
"
);
}
fprintf
(
f
,
"</FieldData>
\n
"
);
fprintf
(
f
,
"</UnstructuredGrid>
\n
"
);
fprintf
(
f
,
"</VTKFile>
\n
"
);
fclose
(
f
);
return
0
;
}
fluid/mesh.h
View file @
eb12c274
...
...
@@ -31,7 +31,7 @@ typedef struct {
int
n_nodes
;
double
*
x
;
int
n_elements
;
u
int
32_t
*
elements
;
int
*
elements
;
int
n_phys
;
char
**
phys_name
;
int
*
phys_n_nodes
;
...
...
python/CMakeLists.txt
View file @
eb12c274
...
...
@@ -26,6 +26,7 @@ set(SRC
lmgc90Interface.py
mesh.py
scontact.py
VTK.py
)
foreach
(
f
${
SRC
}
)
...
...
python/VTK.py
0 → 100644
View file @
eb12c274
import
zlib
import
struct
import
os
import
numpy
as
np
import
re
from
base64
import
b64decode
,
b64encode
from
xml.etree
import
ElementTree
as
ET
def
_typename
(
a
)
:
if
a
.
dtype
==
np
.
int32
:
return
b
"Int32"
elif
a
.
dtype
==
np
.
int64
:
return
b
"Int64"
elif
a
.
dtype
==
np
.
float64
:
return
b
"Float64"
elif
a
.
dtype
==
np
.
float32
:
return
b
"Float32"
else
:
raise
(
NameError
(
"unown array type : "
+
str
(
a
.
dtype
)))
def
_write_array_ascii
(
f
,
a
,
attr
)
:
f
.
write
(
b
'<DataArray %s type="%s" format="ascii">
\n
'
%
(
attr
,
_typename
(
a
)))
f
.
write
((
" "
.
join
(
map
(
str
,
a
.
flat
))).
encode
(
"utf-8"
))
f
.
write
(
b
'
\n
</DataArray>
\n
'
)
def
_write_array
(
f
,
appended
,
a
,
ds
,
attr
)
:
f
.
write
(
b
'<DataArray %s type="%s" format="appended" offset="%i"/>
\n
'
%
(
attr
,
_typename
(
a
),
len
(
appended
)
-
1
))
data
=
zlib
.
compress
(
a
,
2
)
# return appended + struct.pack("iiii",1,a.size*ds,0,len(data)) + data
return
appended
+
b64encode
(
struct
.
pack
(
"iiii"
,
1
,
a
.
size
*
ds
,
0
,
len
(
data
)))
+
b64encode
(
data
)
def
write
(
basename
,
i
,
t
,
elements
,
x
,
data
,
field_data
)
:
appended
=
b
"_"
filename
=
"%s_%05d.vtu"
%
(
basename
,
i
)
f
=
open
(
filename
,
"wb"
)
f
.
write
(
b
'<?xml version="1.0"?>
\n
'
)
f
.
write
(
b
'<VTKFile type="UnstructuredGrid" format="0.1" compressor="vtkZLibDataCompressor">
\n
'
)
f
.
write
(
b
'<UnstructuredGrid>
\n
'
)
nel
=
elements
.
shape
[
0
]
nptbyel
=
elements
.
shape
[
1
]
npt
=
x
.
shape
[
0
]
f
.
write
(
b
'<Piece NumberOfPoints="%i" NumberOfCells="%i">
\n
'
%
(
npt
,
nel
))
f
.
write
(
b
'<Points>
\n
'
)
appended
=
_write_array
(
f
,
appended
,
x
,
8
,
b
'NumberOfComponents="3"'
)
f
.
write
(
b
'</Points>
\n
'
)
f
.
write
(
b
'<Cells>
\n
'
)
appended
=
_write_array
(
f
,
appended
,
elements
,
4
,
b
'Name="connectivity"'
)
offsets
=
np
.
array
(
range
(
nptbyel
,(
nel
+
1
)
*
nptbyel
,
nptbyel
),
np
.
int32
)
appended
=
_write_array
(
f
,
appended
,
offsets
,
4
,
b
'Name="offsets"'
)
types
=
np
.
ones
(
nel
,
np
.
int32
)
*
(
5
if
elements
.
shape
[
1
]
==
3
else
13
)
appended
=
_write_array
(
f
,
appended
,
types
,
4
,
b
'Name="types"'
)
f
.
write
(
b
'</Cells>
\n
'
)
f
.
write
(
b
'<PointData>
\n
'
)
for
name
,
v
in
data
:
vc
=
np
.
ascontiguousarray
(
v
)
appended
=
_write_array
(
f
,
appended
,
vc
,
8
,
b
'Name="%s" NumberOfComponents="%i"'
%
(
name
.
encode
(
"utf8"
),
vc
.
shape
[
-
1
]))
f
.
write
(
b
'</PointData>
\n
'
)
f
.
write
(
b
'</Piece>
\n
'
)
f
.
write
(
b
'<FieldData>
\n
'
);
for
name
,
v
in
field_data
:
vc
=
np
.
ascontiguousarray
(
v
)
appended
=
_write_array
(
f
,
appended
,
vc
,
4
,
b
'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"'
%
(
name
,
vc
.
shape
[
0
],
vc
.
shape
[
1
]))
f
.
write
(
b
'</FieldData>
\n
'
);
f
.
write
(
b
'</UnstructuredGrid>
\n
'
)
f
.
write
(
b
'<AppendedData encoding="base64">
\n
'
)
f
.
write
(
appended
)
f
.
write
(
b
'</AppendedData>
\n
'
)
f
.
write
(
b
'</VTKFile>
\n
'
)
f
.
close
()
def
_get_array_info
(
a
)
:
ncomp
=
int
(
a
[
"NumberOfComponents"
])
if
"NumberOfComponents"
in
a
else
None
ntuples
=
int
(
a
[
"NumberOfTuples"
])
if
"NumberOfTuples"
in
a
else
None
return
(
a
[
"Name"
],
a
[
"type"
],
int
(
a
[
"offset"
]),(
ntuples
,
ncomp
))
def
_get_binary_array
(
appended
,
pos
,
dtype
,
shape
)
:
data
=
appended
[
pos
:
pos
+
24
]
_
,
n
,
_
,
s
=
struct
.
unpack
(
"iiii"
,
b64decode
(
data
))
data
=
appended
[
pos
+
24
:
pos
+
24
-
(
-
s
//
3
)
*
4
]
data
=
b64decode
(
data
)
data
=
zlib
.
decompress
(
data
)
return
np
.
frombuffer
(
data
,
dtype
).
reshape
(
shape
)
def
read
(
fname
)
:
f
=
open
(
fname
,
"rb"
)
root
=
ET
.
parse
(
f
).
getroot
()
grid
=
root
.
find
(
"UnstructuredGrid"
)
piece
=
grid
.
find
(
"Piece"
)
npoints
=
int
(
piece
.
attrib
[
"NumberOfPoints"
])
ncells
=
int
(
piece
.
attrib
[
"NumberOfCells"
])
appended
=
root
.
find
(
"AppendedData"
).
text
shift
=
appended
.
find
(
"_"
)
+
1
dim
=
2
pointsoffset
=
int
(
piece
.
find
(
"Points/DataArray"
).
attrib
[
"offset"
])
x
=
_get_binary_array
(
appended
,
shift
+
pointsoffset
,
np
.
float64
,(
npoints
,
3
))
for
f
in
piece
.
find
(
"Cells"
)
:
if
f
.
attrib
[
"Name"
]
==
"connectivity"
:
connectivityoffset
=
int
(
f
.
attrib
[
"offset"
])
el
=
_get_binary_array
(
appended
,
shift
+
connectivityoffset
,
np
.
int32
,(
ncells
,
dim
+
1
))
fdata
=
[]
for
f
in
grid
.
find
(
"FieldData"
)
:
name
,
ntype
,
offset
,
shape
=
_get_array_info
(
f
.
attrib
)
v
=
_get_binary_array
(
appended
,
shift
+
offset
,
np
.
int32
,
shape
)
fdata
.
append
((
name
,
v
))
data
=
[]
for
f
in
piece
.
find
(
"PointData"
)
:
name
,
ntype
,
offset
,
shape
=
_get_array_info
(
f
.
attrib
)
v
=
_get_binary_array
(
appended
,
shift
+
offset
,
np
.
float64
,(
npoints
,
shape
[
1
]))
data
.
append
((
name
,
v
))
return
x
,
el
,
data
,
fdata
python/fluid.py
View file @
eb12c274
...
...
@@ -32,6 +32,7 @@ from ctypes import *
import
numpy
import
signal
import
os
from
.
import
VTK
dir_path
=
os
.
path
.
dirname
(
os
.
path
.
realpath
(
__file__
))
lib2
=
CDLL
(
os
.
path
.
join
(
dir_path
,
"libmbfluid2.so"
))
...
...
@@ -39,7 +40,6 @@ lib3 = CDLL(os.path.join(dir_path,"libmbfluid3.so"))
signal
.
signal
(
signal
.
SIGINT
,
signal
.
SIG_DFL
)
BNDCB
=
CFUNCTYPE
(
None
,
c_int
,
POINTER
(
c_double
),
POINTER
(
c_double
))
class
FluidProblem
:
...
...
@@ -83,7 +83,7 @@ class FluidProblem :
def
__del__
(
self
):
"""Delete the fluid structure."""
if
(
self
.
_fp
!=
None
)
: