Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
D
dune-elasticity
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Requirements
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Locked files
Build
Pipelines
Jobs
Pipeline schedules
Test cases
Artifacts
Deploy
Releases
Package registry
Container registry
Model registry
Operate
Environments
Terraform modules
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Code review analytics
Issue analytics
Insights
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
This is an archived project. Repository and other project resources are read-only.
Show more breadcrumbs
Ansgar Burchardt
dune-elasticity
Commits
7a2b5152
Commit
7a2b5152
authored
6 years ago
by
Jonathan Youett
Browse files
Options
Downloads
Patches
Plain Diff
Major cleanup
parent
ebaf0285
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
test/adolcmaterialtest.cc
+38
-76
38 additions, 76 deletions
test/adolcmaterialtest.cc
test/materialtest.cc
+152
-164
152 additions, 164 deletions
test/materialtest.cc
with
190 additions
and
240 deletions
test/adolcmaterialtest.cc
+
38
−
76
View file @
7a2b5152
...
@@ -14,25 +14,19 @@
...
@@ -14,25 +14,19 @@
#include
<dune/istl/matrixindexset.hh>
#include
<dune/istl/matrixindexset.hh>
#include
<dune/grid/uggrid.hh>
#include
<dune/grid/uggrid.hh>
#include
<dune/grid/geometrygrid/grid.hh>
#include
<dune/fufem/makesphere.hh>
#include
<dune/fufem/makesphere.hh>
#include
<dune/fufem/makering.hh>
#include
<dune/fufem/makering.hh>
#include
<dune/fufem/functions/basisgridfunction.hh>
#include
<dune/fufem/functions/basisgridfunction.hh>
#include
<dune/fufem/functions/deformationfunction.hh>
#include
<dune/fufem/functionspacebases/p1nodalbasis.hh>
#include
<dune/fufem/functionspacebases/p1nodalbasis.hh>
#include
<dune/fufem/assemblers/operatorassembler.hh>
#include
<dune/fufem/assemblers/operatorassembler.hh>
#include
<dune/fufem/assemblers/functionalassembler.hh>
#include
<dune/fufem/assemblers/functionalassembler.hh>
#include
<dune/fufem/assemblers/localassemblers/adolclinearizationassembler.hh>
#include
<dune/fufem/assemblers/localassemblers/adolchessianassembler.hh>
#define LAURSEN
#define LAURSEN
#include
<dune/elasticity/materials/adolcmaterial.hh>
#include
<dune/elasticity/materials/adolcmaterial.hh>
#include
<dune/elasticity/materials/neohookeanmaterial.hh>
#include
<dune/elasticity/materials/neohookeanmaterial.hh>
#include
<dune/elasticity/materials/adolcneohookeanmaterial.hh>
#include
<dune/elasticity/materials/mooneyrivlinmaterial.hh>
//#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
#include
<dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
#include
<dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
...
@@ -42,116 +36,85 @@ const int dim = 3;
...
@@ -42,116 +36,85 @@ const int dim = 3;
using
namespace
Dune
;
using
namespace
Dune
;
int
main
(
int
argc
,
char
*
argv
[])
try
int
main
(
int
argc
,
char
*
argv
[])
try
{
{
using
FVector
=
FieldVector
<
double
,
dim
>
;
double
E
=
1e5
;
if
(
argc
==
2
)
E
=
atof
(
argv
[
1
]);
typedef
FieldVector
<
double
,
dim
>
FVector
;
typedef
FieldMatrix
<
double
,
dim
,
dim
>
FMatrix
;
typedef
BlockVector
<
FVector
>
VectorType
;
typedef
BCRSMatrix
<
FMatrix
>
MatrixType
;
// ///////////////////////////////////////
// Create the grid
// Create the grid
// ///////////////////////////////////////
using
Grid
=
UGGrid
<
dim
>
;
auto
grid
=
makeSphere
<
Grid
>
(
FVector
(
0
),
0.2
);
typedef
UGGrid
<
dim
>
GridType
;
auto
grid
=
makeSphere
<
GridType
>
(
FVector
(
0
),
0.2
);
grid
->
globalRefine
(
3
);
grid
->
globalRefine
(
3
);
typedef
GridType
::
LeafGridView
GridView
;
GridView
gridView
=
grid
->
leafGridView
();
typedef
P1NodalBasis
<
GridView
>
FEBasis
;
using
FEBasis
=
P1NodalBasis
<
typename
Grid
::
Leaf
GridView
>
;
FEBasis
feBasis
(
gridView
);
FEBasis
feBasis
(
grid
->
leafGrid
View
()
);
// //////////////////////////
// Initial iterate
// Initial iterate
// //////////////////////////
VectorType
x
(
feBasis
.
size
());
double
lower_bound
=
0
;
double
lower_bound
=
0
;
double
upper_bound
=
0.05
;
double
upper_bound
=
0.05
;
std
::
uniform_real_distribution
<
double
>
unif
(
lower_bound
,
upper_bound
);
std
::
uniform_real_distribution
<
double
>
unif
(
lower_bound
,
upper_bound
);
std
::
default_random_engine
re
;
std
::
default_random_engine
re
;
for
(
size_t
i
=
0
;
i
<
x
.
size
();
i
++
)
using
BVector
=
BlockVector
<
FVector
>
;
for
(
int
j
=
0
;
j
<
dim
;
j
++
)
BVector
x
(
feBasis
.
size
());
x
[
i
][
j
]
=
unif
(
re
);
for
(
size_t
i
=
0
;
i
<
x
.
size
();
++
i
)
for
(
int
j
=
0
;
j
<
dim
;
++
j
)
auto
xGridFunc
=
std
::
make_shared
<
BasisGridFunction
<
FEBasis
,
VectorType
>
>
(
feBasis
,
x
);
x
[
i
][
j
]
=
unif
(
re
);
auto
xGridFunc
=
std
::
make_shared
<
BasisGridFunction
<
FEBasis
,
BVector
>
>
(
feBasis
,
x
);
// ////////////////////////////////////////////////////////////
// Create an assembler for the energy functional
// ////////////////////////////////////////////////////////////
double
nu
=
0.4
;
typedef
FEBasis
::
LocalFiniteElement
Lfe
;
// make a material
// make a material
typedef
NeoHookeanMaterial
<
FEBasis
>
MaterialType
;
using
Material
=
NeoHookeanMaterial
<
FEBasis
>
;
Material
Type
material
(
feBasis
,
E
,
nu
);
Material
material
(
feBasis
,
1e5
,
0.4
);
// check against the Adol-C material
// check against the Adol-C material
LocalNeoHookeanEnergy
<
GridType
,
Lfe
>
localNeoHookeEnergy
(
E
,
nu
);
AdolcMaterial
<
FEBasis
>
adolcMaterial
(
feBasis
,
material
,
false
);
typedef
AdolcMaterial
<
FEBasis
>
AdolcMaterialType
;
AdolcMaterialType
adolcMaterial
(
feBasis
,
localNeoHookeEnergy
,
false
);
auto
&
linPaperAssembler
=
material
.
firstDerivative
(
xGridFunc
);
auto
&
linPaperAssembler
=
material
.
firstDerivative
(
xGridFunc
);
auto
&
linAdolcAssembler
=
adolcMaterial
.
firstDerivative
(
xGridFunc
);
auto
&
linAdolcAssembler
=
adolcMaterial
.
firstDerivative
(
xGridFunc
);
FunctionalAssembler
<
FEBasis
>
functionalAssembler
(
feBasis
);
FunctionalAssembler
<
FEBasis
>
functionalAssembler
(
feBasis
);
Vector
Type
adolcGradient
(
feBasis
.
size
());
B
Vector
adolcGradient
(
feBasis
.
size
());
Vector
Type
paperGrad
(
feBasis
.
size
());
B
Vector
paperGrad
(
feBasis
.
size
());
adolcGradient
=
0
;
adolcGradient
=
0
;
paperGrad
=
0
;
paperGrad
=
0
;
Dune
::
Timer
time
;
Timer
time
;
functionalAssembler
.
assemble
(
linAdolcAssembler
,
adolcGradient
);
functionalAssembler
.
assemble
(
linAdolcAssembler
,
adolcGradient
);
std
::
cout
<<
"ADOL-C functional assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
std
::
cout
<<
"ADOL-C functional assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
time
.
reset
();
time
.
reset
();
time
.
start
();
time
.
start
();
functionalAssembler
.
assemble
(
linPaperAssembler
,
paperGrad
);
functionalAssembler
.
assemble
(
linPaperAssembler
,
paperGrad
);
std
::
cout
<<
"Paper&Pen functional assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
std
::
cout
<<
"Paper&Pen functional assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
for
(
size_t
i
=
0
;
i
<
adolcGradient
.
size
();
i
++
)
{
for
(
size_t
i
=
0
;
i
<
adolcGradient
.
size
();
++
i
)
{
FVector
diff
=
adolcGradient
[
i
];
auto
diff
=
adolcGradient
[
i
]
-
paperGrad
[
i
];
diff
-=
paperGrad
[
i
];
if
(
diff
.
two_norm
()
>
1e-
8
)
if
(
diff
.
two_norm
()
>
1e-
6
)
DUNE_THROW
(
Dune
::
Exception
,
"Wrong local derivative, error is "
<<
diff
.
two_norm
());
DUNE_THROW
(
Dune
::
Exception
,
"Wrong local derivative, error is "
<<
diff
.
two_norm
());
}
}
// ////////////////////////
// Test hessian assembler
// ////////////////////////
// Test hessian assembler
auto
&
hessPaperAssembler
=
material
.
secondDerivative
(
xGridFunc
);
auto
&
hessPaperAssembler
=
material
.
secondDerivative
(
xGridFunc
);
auto
&
hessAdolcAssembler
=
adolcMaterial
.
secondDerivative
(
xGridFunc
);
auto
&
hessAdolcAssembler
=
adolcMaterial
.
secondDerivative
(
xGridFunc
);
OperatorAssembler
<
FEBasis
,
FEBasis
>
operatorAssembler
(
feBasis
,
feBasis
);
OperatorAssembler
<
FEBasis
,
FEBasis
>
operatorAssembler
(
feBasis
,
feBasis
);
MatrixType
adolcHessian
,
paperHessian
;
using
Matrix
=
BCRSMatrix
<
FieldMatrix
<
double
,
dim
,
dim
>
>
;
Matrix
adolcHessian
,
paperHessian
;
time
.
reset
();
time
.
reset
();
time
.
start
();
time
.
start
();
operatorAssembler
.
assemble
(
hessAdolcAssembler
,
adolcHessian
);
operatorAssembler
.
assemble
(
hessAdolcAssembler
,
adolcHessian
);
std
::
cout
<<
"ADOL-C operator assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
std
::
cout
<<
"ADOL-C operator assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
time
.
reset
();
time
.
reset
();
time
.
start
();
time
.
start
();
operatorAssembler
.
assemble
(
hessPaperAssembler
,
paperHessian
);
operatorAssembler
.
assemble
(
hessPaperAssembler
,
paperHessian
);
std
::
cout
<<
"Paper&Pen operator assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
std
::
cout
<<
"Paper&Pen operator assembler needed "
<<
time
.
stop
()
<<
std
::
endl
;
for
(
size_t
i
=
0
;
i
<
adolcHessian
.
N
();
i
++
)
{
for
(
size_t
i
=
0
;
i
<
adolcHessian
.
N
();
++
i
)
{
const
auto
&
adolcRow
=
adolcHessian
[
i
];
const
auto
&
adolcRow
=
adolcHessian
[
i
];
const
auto
&
paperRow
=
paperHessian
[
i
];
const
auto
&
paperRow
=
paperHessian
[
i
];
...
@@ -168,19 +131,18 @@ int main (int argc, char *argv[]) try
...
@@ -168,19 +131,18 @@ int main (int argc, char *argv[]) try
DUNE_THROW
(
Dune
::
Exception
,
"Not the same sparsity pattern!"
<<
adolcIt
.
index
()
DUNE_THROW
(
Dune
::
Exception
,
"Not the same sparsity pattern!"
<<
adolcIt
.
index
()
<<
"!="
<<
paperIt
.
index
());
<<
"!="
<<
paperIt
.
index
());
FMatrix
diff
=
*
adolcIt
;
auto
diff
=
*
adolcIt
;
diff
-=
*
paperIt
;
diff
-=
*
paperIt
;
if
(
diff
.
frobenius_norm
()
>
1e-
8
)
if
(
diff
.
frobenius_norm
()
>
1e-
6
)
DUNE_THROW
(
Dune
::
Exception
,
"Wrong local hessian, error is "
<<
diff
.
frobenius_norm
());
DUNE_THROW
(
Dune
::
Exception
,
"Wrong local hessian, error is "
<<
diff
.
frobenius_norm
());
}
}
assert
(
paperIt
==
paperEndIt
);
assert
(
paperIt
==
paperEndIt
);
}
}
// //////////////////////////////
return
0
;
}
catch
(
Exception
e
)
{
}
catch
(
Exception
e
)
{
std
::
cout
<<
e
<<
std
::
endl
;
std
::
cout
<<
e
<<
std
::
endl
;
return
1
;
}
}
This diff is collapsed.
Click to expand it.
test/materialtest.cc
+
152
−
164
View file @
7a2b5152
#include
<config.h>
#include
<config.h>
#include
<iomanip>
#include
<iomanip>
#include
<random>
#include
<random>
#include
<dune/common/fvector.hh>
#include
<dune/fufem/utilities/adolcnamespaceinjections.hh>
#include
<dune/common/fmatrix.hh>
#include
<dune/common/fmatrix.hh>
#include
<dune/common/fvector.hh>
#include
<dune/istl/bvector.hh>
#include
<dune/istl/bcrsmatrix.hh>
#include
<dune/istl/bcrsmatrix.hh>
#include
<dune/istl/bvector.hh>
#include
<dune/istl/matrixindexset.hh>
#include
<dune/istl/matrixindexset.hh>
#include
<dune/grid/uggrid.hh>
#include
<dune/grid/geometrygrid/grid.hh>
#include
<dune/grid/geometrygrid/grid.hh>
#include
<dune/grid/uggrid.hh>
#include
<dune/fufem/assemblers/functionalassembler.hh>
#include
<dune/fufem/makesphere.hh>
#include
<dune/fufem/assemblers/operatorassembler.hh>
#include
<dune/fufem/makering.hh>
#include
<dune/fufem/functions/basisgridfunction.hh>
#include
<dune/fufem/functions/basisgridfunction.hh>
#include
<dune/fufem/functions/deformationfunction.hh>
#include
<dune/fufem/functions/deformationfunction.hh>
#include
<dune/fufem/functionspacebases/p1nodalbasis.hh>
#include
<dune/fufem/functionspacebases/p1nodalbasis.hh>
#include
<dune/fufem/
assemblers/operatorassembler
.hh>
#include
<dune/fufem/
makering
.hh>
#include
<dune/fufem/
assemblers/functionalassembler
.hh>
#include
<dune/fufem/
makesphere
.hh>
#define LAURSEN
#define LAURSEN
#include
<dune/elasticity/materials/neohookeanmaterial.hh>
#include
<dune/elasticity/materials/ogdenmaterial.hh>
#include
<dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
#include
<dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
#include
<dune/elasticity/materials/mooneyrivlinmaterial.hh>
#include
<dune/elasticity/materials/mooneyrivlinmaterial.hh>
#include
<dune/elasticity/materials/neohookeanmaterial.hh>
#include
<dune/elasticity/materials/adolcmaterial.hh>
#include
<dune/elasticity/materials/ogdenmaterial.hh>
template
<
class
Basis
,
class
Material
>
template
<
class
Basis
,
class
Material
>
bool
checkHessian
(
const
Basis
&
basis
,
Material
&
material
,
double
eps
)
bool
checkHessian
(
const
Basis
&
basis
,
Material
&
material
,
double
eps
)
{
{
const
int
dim
=
Basis
::
GridView
::
dimension
;
typedef
Dune
::
BlockVector
<
Dune
::
FieldVector
<
double
,
dim
>
>
VectorType
;
typedef
Dune
::
BCRSMatrix
<
Dune
::
FieldMatrix
<
double
,
dim
,
dim
>
>
MatrixType
;
VectorType
conf
(
basis
.
getGridView
().
size
(
dim
));
conf
=
0
;
double
lower_bound
=
0
;
double
upper_bound
=
0.005
;
std
::
uniform_real_distribution
<
double
>
unif
(
lower_bound
,
upper_bound
);
std
::
default_random_engine
re
;
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
i
++
)
constexpr
int
dim
=
Basis
::
GridView
::
dimension
;
for
(
int
j
=
0
;
j
<
dim
;
j
++
)
using
Vector
=
Dune
::
BlockVector
<
Dune
::
FieldVector
<
double
,
dim
>>
;
conf
[
i
][
j
]
=
unif
(
re
);
double
lower_bound
{
0
},
upper_bound
{
0.005
};
std
::
uniform_real_distribution
<
double
>
unif
(
lower_bound
,
upper_bound
);
std
::
default_random_engine
re
;
typedef
BasisGridFunction
<
Basis
,
VectorType
>
GridFunc
;
Vector
conf
(
basis
.
getGridView
().
size
(
dim
));
auto
confG
=
std
::
make_shared
<
GridFunc
>
(
basis
,
conf
);
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
++
i
)
for
(
int
j
=
0
;
j
<
dim
;
++
j
)
conf
[
i
][
j
]
=
unif
(
re
);
// assemble Hessian using the material assemblers
using
GridFunc
=
BasisGridFunction
<
Basis
,
Vector
>
;
auto
&
locHessian
=
material
.
secondDerivative
(
conf
G
);
auto
confG
=
std
::
make_shared
<
GridFunc
>
(
basis
,
conf
);
MatrixType
exactHessian
;
// assemble Hessian using the material assemblers
OperatorAssembler
<
Basis
,
Basis
>
opAss
(
basis
,
basis
);
auto
&
locHessian
=
material
.
secondDerivative
(
confG
);
opAss
.
assemble
(
locHessian
,
exactHessian
);
// compute finite difference approximation
OperatorAssembler
<
Basis
,
Basis
>
opAss
(
basis
,
basis
);
// dense matrix for the FD approximation
using
Matrix
=
Dune
::
BCRSMatrix
<
Dune
::
FieldMatrix
<
double
,
dim
,
dim
>>
;
Dune
::
MatrixIndexSet
fdIndex
(
conf
.
size
(),
conf
.
size
());
Matrix
exactHessian
;
opAss
.
assemble
(
locHessian
,
exactHessian
);
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
i
++
)
// compute finite difference approximation
for
(
size_t
j
=
0
;
j
<
conf
.
size
();
j
++
)
fdIndex
.
add
(
i
,
j
);
// now the values
// dense matrix for the FD approximation
MatrixType
fdMat
;
Dune
::
MatrixIndexSet
fdIndex
(
conf
.
size
(),
conf
.
size
());
fdIndex
.
exportIdx
(
fdMat
);
fdMat
=
0
;
FunctionalAssembler
<
Basis
>
funcAss
(
basis
);
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
i
++
)
for
(
size_t
j
=
0
;
j
<
conf
.
size
();
++
j
)
fdIndex
.
add
(
i
,
j
);
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
i
++
)
Matrix
fdMat
;
for
(
int
j
=
0
;
j
<
dim
;
j
++
)
{
fdIndex
.
exportIdx
(
fdMat
);
fdMat
=
0
;
FunctionalAssembler
<
Basis
>
funcAss
(
basis
);
VectorType
forward
=
conf
;
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
++
i
)
VectorType
backward
=
conf
;
for
(
int
j
=
0
;
j
<
dim
;
++
j
)
{
forward
[
i
][
j
]
+=
eps
;
auto
forward
=
conf
;
backward
[
i
][
j
]
-=
eps
;
auto
backward
=
conf
;
std
::
shared_ptr
<
GridFunc
>
fdConf
(
new
GridFunc
(
basis
,
forward
));
forward
[
i
][
j
]
+=
eps
;
backward
[
i
][
j
]
-=
eps
;
// assemble Hessian using the material assemblers
auto
fdConf
=
std
::
make_shared
<
GridFunc
>
(
basis
,
forward
);
typename
Material
::
LocalLinearization
&
locLin
=
material
.
firstDerivative
(
fdConf
);
Vector
Type
forw
;
Vector
forw
;
funcAss
.
assemble
(
locLin
,
forw
);
funcAss
.
assemble
(
material
.
firstDerivative
(
fdConf
),
forw
);
std
::
shared_ptr
<
GridFunc
>
fdConfB
(
new
GridFunc
(
basis
,
backward
));
fdConf
=
std
::
make_shared
<
GridFunc
>
(
basis
,
backward
);
auto
&
locLin2
=
material
.
firstDerivative
(
fdConfB
);
Vector
Type
backw
;
Vector
backw
;
funcAss
.
assemble
(
locLin2
,
backw
);
funcAss
.
assemble
(
material
.
firstDerivative
(
fdConf
),
backw
);
for
(
size_t
k
=
0
;
k
<
forw
.
size
();
k
++
)
for
(
size_t
k
=
0
;
k
<
forw
.
size
();
++
k
)
for
(
int
l
=
0
;
l
<
dim
;
l
++
)
for
(
int
l
=
0
;
l
<
dim
;
++
l
)
fdMat
[
k
][
i
][
l
][
j
]
=
(
forw
[
k
][
l
]
-
backw
[
k
][
l
])
/
(
2
*
eps
);
fdMat
[
k
][
i
][
l
][
j
]
=
(
forw
[
k
][
l
]
-
backw
[
k
][
l
])
/
(
2
*
eps
);
}
}
MatrixType
diff
=
fdMat
;
auto
diff
=
fdMat
;
diff
-=
exactHessian
;
diff
-=
exactHessian
;
std
::
cout
<<
"Maximal difference "
<<
std
::
setprecision
(
12
)
<<
diff
.
infinity_norm
()
<<
"
\n
"
;
std
::
cout
<<
"Maximal difference "
<<
std
::
setprecision
(
12
)
<<
diff
.
infinity_norm
()
<<
"
\n
"
;
return
(
diff
.
infinity_norm
()
<
1e-3
);
return
(
diff
.
infinity_norm
()
<
1e-3
);
}
}
template
<
class
Basis
,
class
Material
>
template
<
class
Basis
,
class
Material
>
bool
checkLinearization
(
const
Basis
&
basis
,
Material
&
material
,
double
eps
)
bool
checkLinearization
(
const
Basis
&
basis
,
Material
&
material
,
double
eps
)
{
{
const
int
dim
=
Basis
::
GridView
::
dimension
;
typedef
Dune
::
BlockVector
<
Dune
::
FieldVector
<
double
,
dim
>
>
VectorType
;
VectorType
conf
(
basis
.
getGridView
().
size
(
dim
))
;
constexpr
int
dim
=
Basis
::
GridView
::
dimension
;
conf
=
0
;
using
Vector
=
Dune
::
BlockVector
<
Dune
::
FieldVector
<
double
,
dim
>
>
;
double
lower_bound
=
0
;
double
lower_bound
{
0
},
upper_bound
{
0.05
};
double
upper_bound
=
0.05
;
std
::
uniform_real_distribution
<
double
>
unif
(
lower_bound
,
upper_bound
);
std
::
uniform_real_distribution
<
double
>
unif
(
lower_bound
,
upper_bound
);
std
::
default_random_engine
re
;
std
::
default_random_engine
re
;
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
i
++
)
Vector
conf
(
basis
.
getGridView
().
size
(
dim
));
for
(
int
j
=
0
;
j
<
dim
;
j
++
)
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
++
i
)
conf
[
i
][
j
]
=
unif
(
re
);
for
(
int
j
=
0
;
j
<
dim
;
++
j
)
conf
[
i
][
j
]
=
unif
(
re
);
typedef
BasisGridFunction
<
Basis
,
Vector
Type
>
GridFunc
;
using
GridFunc
=
BasisGridFunction
<
Basis
,
Vector
>
;
auto
confG
=
std
::
make_shared
<
GridFunc
>
(
basis
,
conf
);
auto
confG
=
std
::
make_shared
<
GridFunc
>
(
basis
,
conf
);
// assemble Hessian using the material assemblers
// assemble Hessian using the material assemblers
auto
&
locLinear
=
material
.
firstDerivative
(
confG
);
auto
&
locLinear
=
material
.
firstDerivative
(
confG
);
VectorType
exactLinear
;
FunctionalAssembler
<
Basis
>
funcAss
(
basis
)
;
FunctionalAssembler
<
Basis
>
funcAss
(
basis
)
;
Vector
exactLinear
;
funcAss
.
assemble
(
locLinear
,
exactLinear
);
funcAss
.
assemble
(
locLinear
,
exactLinear
);
// compute finite difference approximation
// compute finite difference approximation
VectorType
fdLin
=
conf
;
auto
fdLin
=
conf
;
VectorType
forward
=
conf
;
auto
forward
=
conf
;
VectorType
backward
=
conf
;
auto
backward
=
conf
;
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
i
++
)
for
(
size_t
i
=
0
;
i
<
conf
.
size
();
++
i
)
for
(
int
j
=
0
;
j
<
dim
;
j
++
)
{
for
(
int
j
=
0
;
j
<
dim
;
++
j
)
{
forward
[
i
][
j
]
+=
eps
;
forward
[
i
][
j
]
+=
eps
;
backward
[
i
][
j
]
-=
eps
;
backward
[
i
][
j
]
-=
eps
;
double
forE
=
material
.
energy
(
std
::
make_shared
<
GridFunc
>
(
basis
,
forward
));
auto
forE
=
material
.
energy
(
std
::
make_shared
<
GridFunc
>
(
basis
,
forward
));
double
backE
=
material
.
energy
(
std
::
make_shared
<
GridFunc
>
(
basis
,
backward
));
auto
backE
=
material
.
energy
(
std
::
make_shared
<
GridFunc
>
(
basis
,
backward
));
fdLin
[
i
][
j
]
=
(
forE
-
backE
)
/
(
2
*
eps
);
fdLin
[
i
][
j
]
=
(
forE
-
backE
)
/
(
2
*
eps
);
forward
[
i
][
j
]
-=
eps
;
forward
[
i
][
j
]
-=
eps
;
backward
[
i
][
j
]
+=
eps
;
backward
[
i
][
j
]
+=
eps
;
}
}
VectorType
diff
=
fdLin
;
auto
diff
=
fdLin
;
diff
-=
exactLinear
;
diff
-=
exactLinear
;
std
::
cout
<<
"Two norm "
<<
std
::
setprecision
(
12
)
<<
diff
.
two_norm
()
<<
"
\n
"
;
std
::
cout
<<
"Two norm "
<<
std
::
setprecision
(
12
)
<<
diff
.
two_norm
()
<<
"
\n
"
;
return
(
diff
.
two_norm
()
<
1e-4
);
return
(
diff
.
two_norm
()
<
1e-4
);
}
}
using
namespace
Dune
;
using
namespace
Dune
;
int
main
(
int
argc
,
char
*
argv
[])
try
int
main
(
int
argc
,
char
*
argv
[])
try
{
{
// parse eps
double
eps
=
1e-6
;
if
(
argc
==
2
)
eps
=
atof
(
argv
[
1
]);
std
::
cout
<<
"eps "
<<
eps
<<
std
::
endl
;
typedef
UGGrid
<
3
>
GridType
;
// parse eps
double
eps
=
1e-6
;
if
(
argc
==
2
)
eps
=
atof
(
argv
[
1
]);
std
::
cout
<<
"eps "
<<
eps
<<
std
::
endl
;
// make sphere grid
using
Grid
=
UGGrid
<
3
>
;
FieldVector
<
double
,
3
>
center
(
0
);
double
radius
=
1
;
auto
grid
=
makeSphere
<
GridType
>
(
center
,
radius
);
// make sphere grid
FieldVector
<
double
,
3
>
center
(
0
);
double
radius
=
1
;
grid
->
setRefinementType
(
GridType
::
COPY
);
auto
grid
=
makeSphere
<
Grid
>
(
center
,
radius
);
for
(
int
i
=
0
;
i
<
1
;
i
++
)
grid
->
globalRefine
(
1
);
// Create the materials
grid
->
setRefinementType
(
Grid
::
COPY
);
typedef
P1NodalBasis
<
GridType
::
LeafGridView
>
P1Basis
;
for
(
int
i
=
0
;
i
<
1
;
i
++
)
P1Basis
p1Basis
(
grid
->
leafGridView
()
);
grid
->
globalRefine
(
1
);
// make a material
// Create the materials
//typedef MooneyRivlinMaterial<P1Basis> MaterialType;
using
P1Basis
=
P1NodalBasis
<
Grid
::
LeafGridView
>
;
typedef
NeoHookeanMaterial
<
P1Basis
>
MaterialType
;
P1Basis
p1Basis
(
grid
->
leafGridView
());
//typedef OgdenMaterial<P1Basis> MaterialType;
//typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
MaterialType
material
;
// make a material
//material.setup(p1Basis, 1e5,0.3, 100);
/*
material
.
setup
(
p1Basis
,
1e5
,
0.3
);
using MoonRivlinMaterial = MooneyRivlinMaterial<P1Basis>;
using MaterialType = AdolcMaterial<P1Basis>;
MoonRivlinMaterial materialM;
materialM.setup(p1Basis, 1e5, 0.3);
MaterialType material(p1Basis, materialM);
*/
bool
passed
=
checkLinearization
(
p1Basis
,
material
,
eps
);
using
MaterialType
=
NeoHookeanMaterial
<
P1Basis
>
;
passed
=
checkHessian
(
p1Basis
,
material
,
eps
)
and
passed
;
// typedef OgdenMaterial<P1Basis> MaterialType;
// typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
std
::
cout
<<
"Checking 2D case!
\n
"
;
MaterialType
material
(
p1Basis
,
1e5
,
0.3
);
//materialM.setup(p1Basis, 1e5, 0.3);
typedef
UGGrid
<
2
>
GridType2
;
bool
passed
=
checkLinearization
(
p1Basis
,
material
,
eps
)
;
FieldVector
<
double
,
2
>
cent2
(
0
)
;
passed
=
checkHessian
(
p1Basis
,
material
,
eps
)
and
passed
;
auto
grid2
=
makeRingSegment2D
<
GridType2
>
(
cent2
,
0.2
,
1
,
0
,
1.5
,
10
,
false
)
;
std
::
cout
<<
"Checking 2D case!
\n
"
;
grid2
->
setRefinementType
(
GridType2
::
COPY
);
using
Grid2
=
UGGrid
<
2
>
;
for
(
int
i
=
0
;
i
<
1
;
i
++
)
FieldVector
<
double
,
2
>
cent2
(
0
);
grid2
->
globalRefine
(
1
);
// Create the materials
auto
grid2
=
makeRingSegment2D
<
Grid2
>
(
cent2
,
0.2
,
1
,
0
,
1.5
,
10
,
false
);
typedef
P1NodalBasis
<
GridType2
::
LeafGridView
>
P1Basis2
;
P1Basis2
p1Basis2
(
grid2
->
leafGridView
());
// make a material
grid2
->
setRefinementType
(
Grid2
::
COPY
);
//typedef MooneyRivlinMaterial<P1Basis2> MaterialType2;
for
(
int
i
=
0
;
i
<
1
;
i
++
)
typedef
NeoHookeanMaterial
<
P1Basis2
>
MaterialType2
;
grid2
->
globalRefine
(
1
);
//typedef OgdenMaterial<P1Basis2> MaterialType2;
//typedef GeomExactStVenantMaterial<P1Basis2> MaterialType2;
MaterialType2
material
2
;
// Create the
material
s
material2
.
setup
(
p1Basis2
,
1e5
,
0.3
)
;
using
P1Basis2
=
P1NodalBasis
<
Grid2
::
LeafGridView
>
;
//material2.setup(p1Basis2, 1e5, 0.3, 100
);
P1Basis2
p1Basis2
(
grid2
->
leafGridView
()
);
passed
=
checkLinearization
(
p1Basis2
,
material2
,
eps
)
and
passed
;
// make a material
passed
=
checkHessian
(
p1Basis2
,
material2
,
eps
)
and
passed
;
//using MaterialType2 = MooneyRivlinMaterial<P1Basis2>;
using
MaterialType2
=
NeoHookeanMaterial
<
P1Basis2
>
;
// typedef OgdenMaterial<P1Basis2> MaterialType2;
// typedef GeomExactStVenantMaterial<P1Basis2> MaterialType2;
return
!
passed
;
//MaterialType2 material2(p1Basis2, 1e5, 0.3);
MaterialType2
material2
(
p1Basis2
,
1e5
,
0.3
);
}
catch
(
Exception
e
)
{
passed
=
checkLinearization
(
p1Basis2
,
material2
,
eps
)
and
passed
;
passed
=
checkHessian
(
p1Basis2
,
material2
,
eps
)
and
passed
;
std
::
cout
<<
e
<<
std
::
endl
;
return
!
passed
;
}
catch
(
Exception
e
)
{
std
::
cout
<<
e
<<
std
::
endl
;
return
1
;
}
}
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