Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
gmsh
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Terms and privacy
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Larry Price
gmsh
Commits
a4e3e1af
Commit
a4e3e1af
authored
14 years ago
by
Axel Modave
Browse files
Options
Downloads
Patches
Plain Diff
pp
parent
f849a9c3
No related branches found
No related tags found
No related merge requests found
Changes
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
Solver/function.cpp
+245
-164
245 additions, 164 deletions
Solver/function.cpp
Solver/function.h
+164
-99
164 additions, 99 deletions
Solver/function.h
with
409 additions
and
263 deletions
Solver/function.cpp
+
245
−
164
View file @
a4e3e1af
...
...
@@ -11,6 +11,7 @@
#include
"Bindings.h"
struct
functionReplaceCache
{
dataCacheMap
*
map
;
std
::
vector
<
dataCacheDouble
*>
toReplace
;
...
...
@@ -18,21 +19,32 @@ struct functionReplaceCache {
};
/*==============================================================================
* class Function (with additionnal class)
*============================================================================*/
// Constructor and destructor
function
::
function
(
int
nbCol
,
bool
invalidatedOnElement
)
:
_nbCol
(
nbCol
),
_invalidatedOnElement
(
invalidatedOnElement
)
{
}
function
::~
function
()
{
}
function
::
function
(
int
nbCol
,
bool
invalidatedOnElement
)
:
_nbCol
(
nbCol
),
_invalidatedOnElement
(
invalidatedOnElement
){}
// Set
void
function
::
addFunctionReplace
(
functionReplace
&
fr
)
{
fr
.
_master
=
this
;
_functionReplaces
.
push_back
(
&
fr
);
}
void
function
::
setArgument
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
)
{
if
(
f
==
NULL
)
if
(
f
==
NULL
)
throw
;
arguments
.
push_back
(
argument
(
v
,
iMap
,
f
));
dependencies
.
insert
(
dependency
(
iMap
,
f
));
for
(
std
::
set
<
dependency
>::
const_iterator
it
=
f
->
dependencies
.
begin
();
it
!=
f
->
dependencies
.
end
();
it
++
)
{
for
(
std
::
set
<
dependency
>::
const_iterator
it
=
f
->
dependencies
.
begin
();
it
!=
f
->
dependencies
.
end
();
it
++
)
{
if
(
it
->
iMap
>
0
&&
iMap
>
0
)
Msg
::
Error
(
"consecutive secondary caches"
);
dependencies
.
insert
(
dependency
(
iMap
+
it
->
iMap
,
it
->
f
));
...
...
@@ -47,6 +59,158 @@ void function::setArgument(fullMatrix<double> &v, const function *f, int iMap) {
}
}
// Time and DT
functionConstant
*
function
::
_timeFunction
=
NULL
;
functionConstant
*
function
::
_dtFunction
=
NULL
;
functionConstant
*
function
::
getTime
()
{
if
(
!
_timeFunction
)
_timeFunction
=
functionConstantNew
(
0.
);
return
_timeFunction
;
}
functionConstant
*
function
::
getDT
()
{
if
(
!
_dtFunction
)
_dtFunction
=
functionConstantNew
(
0.
);
return
_dtFunction
;
}
// Get informations
functionSolution
*
functionSolution
::
_instance
=
NULL
;
function
*
function
::
getSolution
()
{
return
functionSolution
::
get
();
}
//------------------------------------
// Get Solution Gradient + Additionnal class
//------------------------------------
class
functionSolutionGradient
:
public
function
{
static
functionSolutionGradient
*
_instance
;
functionSolutionGradient
()
:
function
(
0
){}
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the gradient of the solution but this algorithm does not provide the gradient of the solution"
);
throw
;
}
static
function
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionSolutionGradient
();
return
_instance
;
}
};
functionSolutionGradient
*
functionSolutionGradient
::
_instance
=
NULL
;
function
*
function
::
getSolutionGradient
()
{
return
functionSolutionGradient
::
get
();
}
//------------------------------------
// Get Parametric Coordinates + Additionnal class
//------------------------------------
class
functionParametricCoordinates
:
public
function
{
static
functionParametricCoordinates
*
_instance
;
functionParametricCoordinates
()
:
function
(
0
,
false
){}
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the parametric coordinates but this algorithm does not provide the parametric coordinates"
);
throw
;
}
static
function
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionParametricCoordinates
();
return
_instance
;
}
};
functionParametricCoordinates
*
functionParametricCoordinates
::
_instance
=
NULL
;
function
*
function
::
getParametricCoordinates
()
{
return
functionParametricCoordinates
::
get
();
}
//------------------------------------
// Get Normals + Additionnal class
//------------------------------------
class
functionNormals
:
public
function
{
static
functionNormals
*
_instance
;
functionNormals
()
:
function
(
0
){}
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the normals coordinates but this algorithm does not provide the normals"
);
throw
;
}
static
function
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionNormals
();
return
_instance
;
}
};
functionNormals
*
functionNormals
::
_instance
=
NULL
;
function
*
function
::
getNormals
()
{
return
functionNormals
::
get
();
}
/*==============================================================================
* class functionReplace
*============================================================================*/
functionReplace
::
functionReplace
()
{
_nChildren
=
0
;
}
void
functionReplace
::
get
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
)
{
bool
allDepFromParent
=
true
;
for
(
std
::
set
<
function
::
dependency
>::
const_iterator
itDep
=
f
->
dependencies
.
begin
();
itDep
!=
f
->
dependencies
.
end
();
itDep
++
){
bool
depFromParent
=
(
_replaced
.
count
(
*
itDep
)
==
0
);
for
(
std
::
set
<
function
::
dependency
>::
const_iterator
itDep2
=
itDep
->
f
->
dependencies
.
begin
();
itDep2
!=
itDep
->
f
->
dependencies
.
end
()
&&
depFromParent
;
itDep2
++
)
depFromParent
&=
(
_replaced
.
count
(
*
itDep2
)
==
0
);
if
(
depFromParent
)
_master
->
dependencies
.
insert
(
*
itDep
);
else
allDepFromParent
=
false
;
}
function
::
dependency
asDep
(
iMap
,
f
);
if
(
allDepFromParent
&&
_replaced
.
count
(
asDep
)
==
0
)
_master
->
dependencies
.
insert
(
asDep
);
_toCompute
.
push_back
(
function
::
argument
(
v
,
iMap
,
f
));
}
void
functionReplace
::
replace
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
)
{
_replaced
.
insert
(
function
::
dependency
(
iMap
,
f
));
_toReplace
.
push_back
(
function
::
argument
(
v
,
iMap
,
f
));
}
void
functionReplace
::
addChild
()
{
_nChildren
++
;
}
void
functionReplace
::
compute
()
{
for
(
unsigned
i
=
0
;
i
<
_toReplace
.
size
();
i
++
)
currentCache
->
toReplace
[
i
]
->
set
();
for
(
unsigned
i
=
0
;
i
<
_toCompute
.
size
();
i
++
)
_toCompute
[
i
].
val
->
setAsProxy
(
currentCache
->
toCompute
[
i
]
->
get
());
};
/*==============================================================================
* class dataCacheDouble
*============================================================================*/
dataCacheDouble
::
dataCacheDouble
(
dataCacheMap
*
m
,
function
*
f
)
:
_cacheMap
(
*
m
),
_value
(
m
->
getNbEvaluationPoints
(),
f
->
getNbCol
()),
_function
(
f
)
{
...
...
@@ -72,8 +236,6 @@ void dataCacheDouble::_eval() {
_valid
=
true
;
}
dataCacheDouble
&
dataCacheMap
::
get
(
const
function
*
f
,
dataCacheDouble
*
caller
)
{
// do I have a cache for this function ?
dataCacheDouble
*&
r
=
_cacheDoubleMap
[
f
];
...
...
@@ -134,6 +296,12 @@ dataCacheDouble &dataCacheMap::get(const function *f, dataCacheDouble *caller) {
return
*
r
;
}
/*==============================================================================
* class dataCacheMap
*============================================================================*/
dataCacheMap
::~
dataCacheMap
()
{
for
(
std
::
set
<
dataCacheDouble
*>::
iterator
it
=
_allDataCaches
.
begin
();
...
...
@@ -145,48 +313,32 @@ dataCacheMap::~dataCacheMap()
}
}
void
functionReplace
::
replace
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
)
{
_replaced
.
insert
(
function
::
dependency
(
iMap
,
f
));
_toReplace
.
push_back
(
function
::
argument
(
v
,
iMap
,
f
));
}
void
functionReplace
::
get
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
)
{
bool
allDepFromParent
=
true
;
for
(
std
::
set
<
function
::
dependency
>::
const_iterator
itDep
=
f
->
dependencies
.
begin
();
itDep
!=
f
->
dependencies
.
end
();
itDep
++
){
bool
depFromParent
=
(
_replaced
.
count
(
*
itDep
)
==
0
);
for
(
std
::
set
<
function
::
dependency
>::
const_iterator
itDep2
=
itDep
->
f
->
dependencies
.
begin
();
itDep2
!=
itDep
->
f
->
dependencies
.
end
()
&&
depFromParent
;
itDep2
++
)
depFromParent
&=
(
_replaced
.
count
(
*
itDep2
)
==
0
);
if
(
depFromParent
)
_master
->
dependencies
.
insert
(
*
itDep
);
else
allDepFromParent
=
false
;
void
dataCacheMap
::
setNbEvaluationPoints
(
int
nbEvaluationPoints
)
{
_nbEvaluationPoints
=
nbEvaluationPoints
;
for
(
std
::
set
<
dataCacheDouble
*>::
iterator
it
=
_allDataCaches
.
begin
();
it
!=
_allDataCaches
.
end
();
it
++
){
(
*
it
)
->
resize
(
nbEvaluationPoints
);
}
function
::
dependency
asDep
(
iMap
,
f
);
if
(
allDepFromParent
&&
_replaced
.
count
(
asDep
)
==
0
)
_master
->
dependencies
.
insert
(
asDep
);
_toCompute
.
push_back
(
function
::
argument
(
v
,
iMap
,
f
));
}
void
functionReplace
::
addChild
()
{
_nChildren
++
;
}
functionReplace
::
functionReplace
(){
_nChildren
=
0
;
for
(
std
::
list
<
dataCacheMap
*>::
iterator
it
=
_children
.
begin
();
it
!=
_children
.
end
();
it
++
)
{
(
*
it
)
->
setNbEvaluationPoints
(
nbEvaluationPoints
);
}
}
void
functionReplace
::
compute
(){
for
(
unsigned
i
=
0
;
i
<
_toReplace
.
size
();
i
++
){
currentCache
->
toReplace
[
i
]
->
set
();
}
for
(
unsigned
i
=
0
;
i
<
_toCompute
.
size
();
i
++
)
{
_toCompute
[
i
].
val
->
setAsProxy
(
currentCache
->
toCompute
[
i
]
->
get
());
}
};
// now some example of functions
// constant values copied over each line
/*==============================================================================
* Some examples of functions
*============================================================================*/
// functionConstant (constant values copied over each line)
void
functionConstant
::
set
(
double
val
)
{
if
(
getNbCol
()
!=
1
)
Msg
::
Error
(
"set scalar value on a vectorial constant function"
);
_source
(
0
,
0
)
=
val
;
}
functionConstant
*
functionConstantNew
(
double
v
)
{
std
::
vector
<
double
>
vec
(
1
);
...
...
@@ -198,25 +350,18 @@ functionConstant *functionConstantNew(const std::vector<double> &v) {
return
new
functionConstant
(
v
);
}
void
functionConstant
::
set
(
double
val
)
{
if
(
getNbCol
()
!=
1
)
{
Msg
::
Error
(
"set scalar value on a vectorial constant function"
);
}
_source
(
0
,
0
)
=
val
;
}
// functionSum
class
functionSum
:
public
function
{
public:
public:
fullMatrix
<
double
>
_f0
,
_f1
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
val
.
size2
();
j
++
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
val
.
size2
();
j
++
)
val
(
i
,
j
)
=
_f0
(
i
,
j
)
+
_f1
(
i
,
j
);
}
}
functionSum
(
const
function
*
f0
,
const
function
*
f1
)
:
function
(
f0
->
getNbCol
()){
functionSum
(
const
function
*
f0
,
const
function
*
f1
)
:
function
(
f0
->
getNbCol
())
{
if
(
f0
->
getNbCol
()
!=
f1
->
getNbCol
())
{
Msg
::
Error
(
"trying to sum 2 functions of different sizes
\n
"
);
throw
;
...
...
@@ -230,16 +375,18 @@ function *functionSumNew(const function *f0, const function *f1) {
return
new
functionSum
(
f0
,
f1
);
}
// functionProd
class
functionProd
:
public
function
{
public:
public:
fullMatrix
<
double
>
_f0
,
_f1
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
val
.
size2
();
j
++
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
val
.
size2
();
j
++
)
val
(
i
,
j
)
=
_f0
(
i
,
j
)
*
_f1
(
i
,
j
);
}
}
functionProd
(
const
function
*
f0
,
const
function
*
f1
)
:
function
(
f0
->
getNbCol
()){
functionProd
(
const
function
*
f0
,
const
function
*
f1
)
:
function
(
f0
->
getNbCol
())
{
if
(
f0
->
getNbCol
()
!=
f1
->
getNbCol
())
{
Msg
::
Error
(
"trying to compute product of 2 functions of different sizes
\n
"
);
throw
;
...
...
@@ -253,15 +400,18 @@ function *functionProdNew(const function *f0, const function *f1) {
return
new
functionProd
(
f0
,
f1
);
}
// functionExtractComp
class
functionExtractComp
:
public
function
{
public:
fullMatrix
<
double
>
_f0
;
double
_iComp
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
val
(
i
,
0
)
=
_f0
(
i
,
_iComp
);
}
functionExtractComp
(
const
function
*
f0
,
const
int
iComp
)
:
function
(
1
){
functionExtractComp
(
const
function
*
f0
,
const
int
iComp
)
:
function
(
1
)
{
setArgument
(
_f0
,
f0
);
_iComp
=
iComp
;
}
...
...
@@ -272,17 +422,18 @@ function *functionExtractCompNew(const function *f0, const int iComp) {
}
// functionScale
class
functionScale
:
public
function
{
public:
public:
fullMatrix
<
double
>
_f0
;
double
_s
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
val
.
size2
();
j
++
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
val
.
size2
();
j
++
)
val
(
i
,
j
)
=
_f0
(
i
,
j
)
*
_s
;
}
}
functionScale
(
const
function
*
f0
,
const
double
s
)
:
function
(
f0
->
getNbCol
()){
functionScale
(
const
function
*
f0
,
const
double
s
)
:
function
(
f0
->
getNbCol
())
{
setArgument
(
_f0
,
f0
);
_s
=
s
;
}
...
...
@@ -292,12 +443,14 @@ function *functionScaleNew(const function *f0, const double s) {
return
new
functionScale
(
f0
,
s
);
}
// get XYZ coordinates
// functionCoordinates (get XYZ coordinates)
class
functionCoordinates
:
public
function
{
static
functionCoordinates
*
_instance
;
fullMatrix
<
double
>
uvw
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
xyz
){
for
(
int
i
=
0
;
i
<
uvw
.
size1
();
i
++
){
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
xyz
)
{
for
(
int
i
=
0
;
i
<
uvw
.
size1
();
i
++
)
{
SPoint3
p
;
m
->
getElement
()
->
pnt
(
uvw
(
i
,
0
),
uvw
(
i
,
1
),
uvw
(
i
,
2
),
p
);
xyz
(
i
,
0
)
=
p
.
x
();
...
...
@@ -305,10 +458,9 @@ class functionCoordinates : public function {
xyz
(
i
,
2
)
=
p
.
z
();
}
}
functionCoordinates
()
:
function
(
3
)
{
functionCoordinates
()
:
function
(
3
)
{
// constructor is private only 1 instance can exists, call get to access the instance
setArgument
(
uvw
,
function
::
getParametricCoordinates
());
};
// constructor is private only 1 instance can exists, call get to access the instance
};
public
:
static
function
*
get
()
{
if
(
!
_instance
)
...
...
@@ -316,81 +468,27 @@ class functionCoordinates : public function {
return
_instance
;
}
};
functionCoordinates
*
functionCoordinates
::
_instance
=
NULL
;
functionSolution
*
functionSolution
::
_instance
=
NULL
;
function
*
function
::
getSolution
()
{
return
functionSolution
::
get
();
}
functionCoordinates
*
functionCoordinates
::
_instance
=
NULL
;
class
functionSolutionGradient
:
public
function
{
static
functionSolutionGradient
*
_instance
;
functionSolutionGradient
()
:
function
(
0
){}
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the gradient of the solution but this algorithm does not provide the gradient of the solution"
);
throw
;
}
static
function
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionSolutionGradient
();
return
_instance
;
}
};
functionSolutionGradient
*
functionSolutionGradient
::
_instance
=
NULL
;
function
*
function
::
getSolutionGradient
()
{
return
functionSolutionGradient
::
get
();
function
*
getFunctionCoordinates
()
{
return
functionCoordinates
::
get
();
}
class
functionParametricCoordinates
:
public
function
{
static
functionParametricCoordinates
*
_instance
;
functionParametricCoordinates
()
:
function
(
0
,
false
){}
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the parametric coordinates but this algorithm does not provide the parametric coordinates"
);
throw
;
}
static
function
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionParametricCoordinates
();
return
_instance
;
}
};
functionParametricCoordinates
*
functionParametricCoordinates
::
_instance
=
NULL
;
function
*
function
::
getParametricCoordinates
()
{
return
functionParametricCoordinates
::
get
();
}
class
functionNormals
:
public
function
{
static
functionNormals
*
_instance
;
functionNormals
()
:
function
(
0
){}
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the normals coordinates but this algorithm does not provide the normals"
);
throw
;
}
static
function
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionNormals
();
return
_instance
;
}
};
functionNormals
*
functionNormals
::
_instance
=
NULL
;
function
*
function
::
getNormals
()
{
return
functionNormals
::
get
();
}
// functionStructuredGridFile
class
functionStructuredGridFile
:
public
function
{
fullMatrix
<
double
>
coord
;
public:
public:
int
n
[
3
];
double
d
[
3
],
o
[
3
];
double
get
(
int
i
,
int
j
,
int
k
){
double
get
(
int
i
,
int
j
,
int
k
)
{
return
v
[(
i
*
n
[
1
]
+
j
)
*
n
[
2
]
+
k
];
}
double
*
v
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
){
for
(
int
pt
=
0
;
pt
<
val
.
size1
();
pt
++
){
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
)
{
for
(
int
pt
=
0
;
pt
<
val
.
size1
();
pt
++
)
{
double
xi
[
3
];
int
id
[
3
];
for
(
int
i
=
0
;
i
<
3
;
i
++
){
...
...
@@ -401,7 +499,7 @@ class functionStructuredGridFile : public function {
xi
[
i
]
=
std
::
min
(
1.
,
std
::
max
(
0.
,
xi
[
i
]));
}
val
(
pt
,
0
)
=
get
(
id
[
0
]
,
id
[
1
]
,
id
[
2
]
)
*
(
1
-
xi
[
0
])
*
(
1
-
xi
[
1
])
*
(
1
-
xi
[
2
])
get
(
id
[
0
]
,
id
[
1
]
,
id
[
2
]
)
*
(
1
-
xi
[
0
])
*
(
1
-
xi
[
1
])
*
(
1
-
xi
[
2
])
+
get
(
id
[
0
]
,
id
[
1
]
,
id
[
2
]
+
1
)
*
(
1
-
xi
[
0
])
*
(
1
-
xi
[
1
])
*
(
xi
[
2
])
+
get
(
id
[
0
]
,
id
[
1
]
+
1
,
id
[
2
]
)
*
(
1
-
xi
[
0
])
*
(
xi
[
1
])
*
(
1
-
xi
[
2
])
+
get
(
id
[
0
]
,
id
[
1
]
+
1
,
id
[
2
]
+
1
)
*
(
1
-
xi
[
0
])
*
(
xi
[
1
])
*
(
xi
[
2
])
...
...
@@ -420,9 +518,8 @@ class functionStructuredGridFile : public function {
input
>>
o
[
0
]
>>
o
[
1
]
>>
o
[
2
]
>>
d
[
0
]
>>
d
[
1
]
>>
d
[
2
]
>>
n
[
0
]
>>
n
[
1
]
>>
n
[
2
];
int
nt
=
n
[
0
]
*
n
[
1
]
*
n
[
2
];
v
=
new
double
[
nt
];
for
(
int
i
=
0
;
i
<
nt
;
i
++
)
{
for
(
int
i
=
0
;
i
<
nt
;
i
++
)
input
>>
v
[
i
];
}
}
else
{
input
.
read
((
char
*
)
o
,
3
*
sizeof
(
double
));
input
.
read
((
char
*
)
d
,
3
*
sizeof
(
double
));
...
...
@@ -432,11 +529,14 @@ class functionStructuredGridFile : public function {
input
.
read
((
char
*
)
v
,
nt
*
sizeof
(
double
));
}
}
~
functionStructuredGridFile
(){
~
functionStructuredGridFile
()
{
delete
[]
v
;
}
};
// functionLua
#ifdef HAVE_LUA
class
functionLua
:
public
function
{
lua_State
*
_L
;
...
...
@@ -454,25 +554,15 @@ class functionLua : public function {
:
function
(
nbCol
),
_luaFunctionName
(
luaFunctionName
),
_L
(
L
)
{
args
.
resize
(
dependencies
.
size
());
for
(
int
i
=
0
;
i
<
dependencies
.
size
();
i
++
)
{
setArgument
(
args
[
i
],
dependencies
[
i
]);
}
for
(
int
i
=
0
;
i
<
dependencies
.
size
();
i
++
)
setArgument
(
args
[
i
],
dependencies
[
i
]);
}
};
#endif
// functionC
void
dataCacheMap
::
setNbEvaluationPoints
(
int
nbEvaluationPoints
)
{
_nbEvaluationPoints
=
nbEvaluationPoints
;
for
(
std
::
set
<
dataCacheDouble
*>::
iterator
it
=
_allDataCaches
.
begin
();
it
!=
_allDataCaches
.
end
();
it
++
){
(
*
it
)
->
resize
(
nbEvaluationPoints
);
}
for
(
std
::
list
<
dataCacheMap
*>::
iterator
it
=
_children
.
begin
();
it
!=
_children
.
end
();
it
++
)
{
(
*
it
)
->
setNbEvaluationPoints
(
nbEvaluationPoints
);
}
}
//functionC
class
functionC
:
public
function
{
std
::
vector
<
fullMatrix
<
double
>
>
args
;
void
(
*
callback
)(
void
);
...
...
@@ -551,6 +641,12 @@ class functionC : public function {
}
};
/*==============================================================================
* Bindings
*============================================================================*/
void
function
::
registerBindings
(
binding
*
b
){
classBinding
*
cb
=
b
->
addClass
<
function
>
(
"Function"
);
cb
->
setDescription
(
"A generic function that can be evaluated on a set of points. Functions can call other functions and their values are cached so that if two different functions call the same function f, f is only evaluated once."
);
...
...
@@ -619,18 +715,3 @@ void function::registerBindings(binding *b){
#endif
}
functionConstant
*
function
::
_timeFunction
=
NULL
;
functionConstant
*
function
::
getTime
()
{
if
(
!
_timeFunction
)
_timeFunction
=
functionConstantNew
(
0.
);
return
_timeFunction
;
}
functionConstant
*
function
::
_dtFunction
=
NULL
;
functionConstant
*
function
::
getDT
()
{
if
(
!
_dtFunction
)
_dtFunction
=
functionConstantNew
(
0.
);
return
_dtFunction
;
}
function
*
getFunctionCoordinates
(){
return
functionCoordinates
::
get
();
}
This diff is collapsed.
Click to expand it.
Solver/function.h
+
164
−
99
View file @
a4e3e1af
...
...
@@ -7,71 +7,162 @@
#include
<list>
#include
<string>
#include
<vector>
class
dataCacheMap
;
class
MElement
;
class
binding
;
class
function
;
class
dataCacheDouble
;
class
dataCacheMap
;
class
dgDataCacheMap
;
class
function
;
class
functionConstant
;
class
functionReplace
;
struct
functionReplaceCache
;
class
MElement
;
class
binding
;
/*==============================================================================
* class function : An abstract interface to functions
*============================================================================*/
// An abstract interface to functions
// more explanation at the head of this file
class
function
{
public
:
int
_nbCol
;
bool
_invalidatedOnElement
;
std
::
vector
<
functionReplace
*>
_functionReplaces
;
public:
// Additionnal types
class
dependency
{
public
:
unsigned
iMap
;
const
function
*
f
;
dependency
(
int
iMap_
,
const
function
*
f_
){
iMap
=
iMap_
;
f
=
f_
;
}
public:
unsigned
iMap
;
const
function
*
f
;
dependency
(
int
iMap_
,
const
function
*
f_
)
{
iMap
=
iMap_
;
f
=
f_
;
}
bool
operator
<
(
const
dependency
&
d
)
const
{
return
(
d
.
iMap
<
iMap
||
d
.
f
<
f
);
}
};
void
printDep
()
{
for
(
std
::
set
<
dependency
>::
iterator
it
=
dependencies
.
begin
();
it
!=
dependencies
.
end
();
it
++
)
printf
(
"%i %p
\n
"
,
it
->
iMap
,
it
->
f
);
}
class
argument
{
//iMap is the id of the dataCacheMap, e.g. on interfaces
public:
int
iMap
;
const
function
*
f
;
fullMatrix
<
double
>
*
val
;
argument
(
fullMatrix
<
double
>
&
v
,
int
iMap_
,
const
function
*
f_
){
val
=
&
v
;
iMap
=
iMap_
;
f
=
f_
;
}
public:
int
iMap
;
// id of the dataCacheMap, e.g. on interfaces
const
function
*
f
;
fullMatrix
<
double
>
*
val
;
argument
(
fullMatrix
<
double
>
&
v
,
int
iMap_
,
const
function
*
f_
)
{
val
=
&
v
;
iMap
=
iMap_
;
f
=
f_
;
}
};
// Data
int
_nbCol
;
bool
_invalidatedOnElement
;
std
::
vector
<
functionReplace
*>
_functionReplaces
;
std
::
vector
<
argument
>
arguments
;
std
::
set
<
dependency
>
dependencies
;
void
setArgument
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
=
0
);
void
addFunctionReplace
(
functionReplace
&
fr
);
private
:
static
functionConstant
*
_timeFunction
;
static
functionConstant
*
_dtFunction
;
public
:
function
(
int
nbCol
,
bool
invalidatedOnElement
=
true
);
virtual
~
function
();
virtual
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
res
)
=
0
;
void
addFunctionReplace
(
functionReplace
&
fr
);
void
setArgument
(
fullMatrix
<
double
>
&
v
,
const
function
*
f
,
int
iMap
=
0
);
virtual
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
res
)
=
0
;
virtual
void
registerInDataCacheMap
(
dataCacheMap
*
m
,
dataCacheDouble
*
d
)
{}
inline
bool
isInvalitedOnElement
()
{
return
_invalidatedOnElement
;}
inline
int
getNbCol
()
const
{
return
_nbCol
;}
// Get or print information
static
void
registerBindings
(
binding
*
b
);
inline
bool
isInvalitedOnElement
()
{
return
_invalidatedOnElement
;}
inline
int
getNbCol
()
const
{
return
_nbCol
;}
static
functionConstant
*
getTime
();
static
functionConstant
*
getDT
();
private
:
static
functionConstant
*
_timeFunction
;
static
functionConstant
*
_dtFunction
;
public
:
static
function
*
getSolution
();
static
function
*
getSolutionGradient
();
static
function
*
getParametricCoordinates
();
static
function
*
getNormals
();
static
functionConstant
*
getTime
();
static
functionConstant
*
getDT
();
void
printDep
()
{
for
(
std
::
set
<
dependency
>::
iterator
it
=
dependencies
.
begin
();
it
!=
dependencies
.
end
();
it
++
)
printf
(
"%i %p
\n
"
,
it
->
iMap
,
it
->
f
);
}
// Bindings
static
void
registerBindings
(
binding
*
b
);
};
//------------------------------------
// Additionnal class to get solution
//------------------------------------
class
functionSolution
:
public
function
{
static
functionSolution
*
_instance
;
functionSolution
()
:
function
(
0
)
{};
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the solution but this algorithm does not provide the solution"
);
throw
;
}
void
setNbFields
(
int
nbFields
)
{
_nbCol
=
nbFields
;
}
static
functionSolution
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionSolution
();
return
_instance
;
}
};
/*==============================================================================
* class functionReplace
*============================================================================*/
class
functionReplace
{
friend
class
dataCacheMap
;
friend
class
dataCacheDouble
;
public:
int
_nChildren
;
function
*
_master
;
functionReplaceCache
*
currentCache
;
std
::
set
<
function
::
dependency
>
_replaced
;
std
::
set
<
function
::
dependency
>
_fromParent
;
std
::
vector
<
function
::
argument
>
_toReplace
;
std
::
vector
<
function
::
argument
>
_toCompute
;
functionReplace
();
void
get
(
fullMatrix
<
double
>
&
v
,
const
function
*
,
int
iMap
=
0
);
void
replace
(
fullMatrix
<
double
>
&
v
,
const
function
*
,
int
iMap
=
0
);
void
compute
();
void
addChild
();
};
// dataCache when the value is a matrix of double
// the user should provide the number of rows by evaluating points and the number of columns
// then the size of the matrix is automatically adjusted
/*==============================================================================
* class dataCacheDouble :
* dataCache when the value is a matrix of double
* the user should provide the number of rows by evaluating points and the number of columns
* then the size of the matrix is automatically adjusted
*============================================================================*/
class
dataCacheDouble
{
friend
class
dataCacheMap
;
friend
class
dgDataCacheMap
;
...
...
@@ -80,14 +171,13 @@ class dataCacheDouble {
function
*
_function
;
dataCacheMap
&
_cacheMap
;
std
::
vector
<
functionReplaceCache
>
functionReplaceCaches
;
protected:
protected:
// pointers to all of the dataCache depending on me
std
::
set
<
dataCacheDouble
*>
_dependOnMe
;
std
::
set
<
dataCacheDouble
*>
_iDependOn
;
fullMatrix
<
double
>
_value
;
bool
_valid
;
dataCacheDouble
(
dataCacheMap
*
,
function
*
f
);
// do the actual computation and put the result into _value
void
_eval
();
void
resize
(
int
nrow
);
...
...
@@ -123,8 +213,11 @@ class dataCacheDouble {
};
class
dgDataCacheMap
;
// more explanation at the head of this file
/*==============================================================================
* class dataCacheMap
*============================================================================*/
class
dataCacheMap
{
public:
dataCacheMap
*
_parent
;
...
...
@@ -134,10 +227,13 @@ class dataCacheMap {
std
::
map
<
const
function
*
,
dataCacheDouble
*>
_cacheDoubleMap
;
std
::
set
<
dataCacheDouble
*>
_allDataCaches
;
std
::
set
<
dataCacheDouble
*>
_toInvalidateOnElement
;
MElement
*
_element
;
void
addDataCacheDouble
(
dataCacheDouble
*
data
,
bool
invalidatedOnElement
){
dataCacheMap
()
{
_nbEvaluationPoints
=
0
;
_parent
=
NULL
;
}
~
dataCacheMap
();
void
addDataCacheDouble
(
dataCacheDouble
*
data
,
bool
invalidatedOnElement
)
{
_allDataCaches
.
insert
(
data
);
if
(
invalidatedOnElement
)
_toInvalidateOnElement
.
insert
(
data
);
...
...
@@ -157,18 +253,14 @@ class dataCacheMap {
dataCacheDouble
&
get
(
const
function
*
f
,
dataCacheDouble
*
caller
=
0
);
virtual
void
setElement
(
MElement
*
element
)
{
_element
=
element
;
for
(
std
::
set
<
dataCacheDouble
*>::
iterator
it
=
_toInvalidateOnElement
.
begin
();
it
!=
_toInvalidateOnElement
.
end
();
it
++
)
{
for
(
std
::
set
<
dataCacheDouble
*>::
iterator
it
=
_toInvalidateOnElement
.
begin
();
it
!=
_toInvalidateOnElement
.
end
();
it
++
)
{
(
*
it
)
->
_valid
=
false
;
}
for
(
std
::
list
<
dataCacheMap
*>::
iterator
it
=
_children
.
begin
();
it
!=
_children
.
end
();
it
++
)
{
for
(
std
::
list
<
dataCacheMap
*>::
iterator
it
=
_children
.
begin
();
it
!=
_children
.
end
();
it
++
)
{
(
*
it
)
->
setElement
(
element
);
}
}
inline
MElement
*
getElement
()
{
return
_element
;}
dataCacheMap
()
{
_nbEvaluationPoints
=
0
;
_parent
=
NULL
;
}
virtual
dataCacheMap
*
newChild
()
{
dataCacheMap
*
m
=
new
dataCacheMap
();
m
->
_parent
=
this
;
...
...
@@ -177,70 +269,43 @@ class dataCacheMap {
return
m
;
}
void
setNbEvaluationPoints
(
int
nbEvaluationPoints
);
inline
int
getNbEvaluationPoints
(){
return
_nbEvaluationPoints
;}
~
dataCacheMap
();
};
class
functionReplace
{
friend
class
dataCacheMap
;
friend
class
dataCacheDouble
;
public
:
function
*
_master
;
int
_nChildren
;
functionReplaceCache
*
currentCache
;
std
::
set
<
function
::
dependency
>
_replaced
;
std
::
set
<
function
::
dependency
>
_fromParent
;
std
::
vector
<
function
::
argument
>
_toReplace
;
std
::
vector
<
function
::
argument
>
_toCompute
;
void
get
(
fullMatrix
<
double
>
&
v
,
const
function
*
,
int
iMap
=
0
);
void
replace
(
fullMatrix
<
double
>
&
v
,
const
function
*
,
int
iMap
=
0
);
void
compute
();
functionReplace
();
void
addChild
();
inline
int
getNbEvaluationPoints
()
{
return
_nbEvaluationPoints
;}
};
functionConstant
*
functionConstantNew
(
const
std
::
vector
<
double
>&
);
functionConstant
*
functionConstantNew
(
double
);
function
*
functionSumNew
(
const
function
*
f0
,
const
function
*
f1
);
function
*
functionProdNew
(
const
function
*
f0
,
const
function
*
f1
);
function
*
functionScaleNew
(
const
function
*
f0
,
const
double
s
);
function
*
functionExtractCompNew
(
const
function
*
f0
,
const
int
iComp
);
class
functionSolution
:
public
function
{
static
functionSolution
*
_instance
;
functionSolution
()
:
function
(
0
){};
// constructor is private only 1 instance can exists, call get to access the instance
public
:
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
sol
)
{
Msg
::
Error
(
"a function requires the solution but this algorithm does not provide the solution"
);
throw
;
}
void
setNbFields
(
int
nbFields
)
{
_nbCol
=
nbFields
;
}
static
functionSolution
*
get
()
{
if
(
!
_instance
)
_instance
=
new
functionSolution
();
return
_instance
;
}
};
/*==============================================================================
* Some example of functions
*============================================================================*/
// functionConstant (constant values copied over each line)
class
functionConstant
:
public
function
{
public:
public:
fullMatrix
<
double
>
_source
;
void
call
(
dataCacheMap
*
m
,
fullMatrix
<
double
>
&
val
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
_source
.
size1
();
j
++
)
{
for
(
int
i
=
0
;
i
<
val
.
size1
();
i
++
)
for
(
int
j
=
0
;
j
<
_source
.
size1
();
j
++
)
val
(
i
,
j
)
=
_source
(
j
,
0
);
}
}
functionConstant
(
std
::
vector
<
double
>
source
)
:
function
(
source
.
size
()){
functionConstant
(
std
::
vector
<
double
>
source
)
:
function
(
source
.
size
())
{
_source
=
fullMatrix
<
double
>
(
source
.
size
(),
1
);
for
(
size_t
i
=
0
;
i
<
source
.
size
();
i
++
)
{
for
(
size_t
i
=
0
;
i
<
source
.
size
();
i
++
)
_source
(
i
,
0
)
=
source
[
i
];
}
}
void
set
(
double
val
);
};
functionConstant
*
functionConstantNew
(
const
std
::
vector
<
double
>&
);
functionConstant
*
functionConstantNew
(
double
);
function
*
functionSumNew
(
const
function
*
f0
,
const
function
*
f1
);
function
*
functionProdNew
(
const
function
*
f0
,
const
function
*
f1
);
function
*
functionScaleNew
(
const
function
*
f0
,
const
double
s
);
function
*
functionExtractCompNew
(
const
function
*
f0
,
const
int
iComp
);
function
*
getFunctionCoordinates
();
#endif
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment