Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
Pries, Jason
Oersted
Commits
01d8f973
Commit
01d8f973
authored
Jul 27, 2017
by
Pries, Jason
Browse files
Finishes removal of Physics templates that are not immediately necessary
FiniteElementMesh remains templated to facilitate storage of variable order elements
parent
38350d71
Changes
10
Hide whitespace changes
Inline
Side-by-side
src/Physics/src/BoundaryCondition.h
View file @
01d8f973
...
...
@@ -17,23 +17,15 @@
class
BoundaryCondition
{
public:
virtual
~
BoundaryCondition
()
{
std
::
cout
<<
"Get rid of unnecessary templates"
<<
std
::
endl
;
};
virtual
void
apply
(
Eigen
::
SparseMatrix
<
double
>
&
bc_matrix
)
const
=
0
;
virtual
void
reduce
(
std
::
set
<
size_t
,
std
::
less
<
size_t
>>
&
index
)
const
=
0
;
};
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
ZeroDirichlet
:
public
BoundaryCondition
{
};
template
<
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
ZeroDirichlet
<
2
,
ElementOrder
,
QuadratureOrder
>
:
public
BoundaryCondition
{
public:
ZeroDirichlet
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
boundaries
)
:
Boundaries
{
boundaries
}
{};
ZeroDirichlet
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
boundaries
)
:
Boundaries
{
boundaries
}
{};
//ZeroDirichlet(std::vector<std::shared_ptr<Curve const>> const &boundaries, FiniteElementMesh<ElementOrder, QuadratureOrder> const &fem) {
ZeroDirichlet
(
std
::
vector
<
std
::
shared_ptr
<
Curve
const
>>
const
&
boundaries
,
std
::
shared_ptr
<
FiniteElementMeshInterface
>
fem
)
{
for
(
auto
const
&
b
:
boundaries
)
{
auto
fb
=
fem
->
boundary
(
b
);
...
...
@@ -52,24 +44,18 @@ public:
};
protected:
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
Boundaries
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
Boundaries
;
};
// TODO: Unify interface between PeriodicBoundary and SlidingBoundary (inherit a common interface/data members)
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
PeriodicBoundaryCondition
:
public
BoundaryCondition
{
};
template
<
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
PeriodicBoundaryCondition
<
2
,
ElementOrder
,
QuadratureOrder
>
:
public
BoundaryCondition
{
public:
PeriodicBoundaryCondition
(
std
::
vector
<
VariableMap
>
map
,
bool
antiperiodic
)
:
Map
{
map
},
Antiperiodic
{
antiperiodic
}
{};
//PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<ElementOrder, QuadratureOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
PeriodicBoundaryCondition
(
std
::
vector
<
PeriodicBoundaryPair
>
const
&
periodic_boundary_pairs
,
std
::
shared_ptr
<
FiniteElementMeshInterface
>
fem
,
bool
antiperiodic
)
:
Antiperiodic
{
antiperiodic
}
{
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b0
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b1
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b0
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b1
;
std
::
vector
<
bool
>
orientation
;
for
(
auto
const
&
pbp
:
periodic_boundary_pairs
)
{
...
...
@@ -81,7 +67,7 @@ public:
init
(
b0
,
b1
,
orientation
);
}
PeriodicBoundaryCondition
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b0
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b1
,
std
::
vector
<
bool
>
orientation
,
bool
antiperiodic
)
:
Antiperiodic
{
antiperiodic
}
{
PeriodicBoundaryCondition
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b0
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b1
,
std
::
vector
<
bool
>
orientation
,
bool
antiperiodic
)
:
Antiperiodic
{
antiperiodic
}
{
init
(
b0
,
b1
,
orientation
);
}
...
...
@@ -124,7 +110,7 @@ protected:
bool
Antiperiodic
;
private:
void
init
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b0
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b1
,
std
::
vector
<
bool
>
orientation
)
{
void
init
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b0
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b1
,
std
::
vector
<
bool
>
orientation
)
{
double_t
sgn
=
sign
();
for
(
size_t
i
=
0
;
i
!=
b0
.
size
();
++
i
)
{
...
...
@@ -156,14 +142,9 @@ private:
}
};
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
SlidingInterface
:
public
BoundaryCondition
{
};
template
<
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
SlidingInterface
<
2
,
ElementOrder
,
QuadratureOrder
>
:
public
BoundaryCondition
{
public:
SlidingInterface
(
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
b0
,
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
b1
,
bool
antiperiodic
=
false
)
:
Antiperiodic
{
antiperiodic
}
{
SlidingInterface
(
std
::
shared_ptr
<
DiscreteBoundary
>
b0
,
std
::
shared_ptr
<
DiscreteBoundary
>
b1
,
bool
antiperiodic
=
false
)
:
Antiperiodic
{
antiperiodic
}
{
std
::
vector
<
size_t
>
const
&
first
=
b0
->
nodes
();
First
=
std
::
vector
<
size_t
>
(
first
.
begin
(),
first
.
end
()
-
1
);
...
...
src/Physics/src/DiscreteBoundary.h
View file @
01d8f973
...
...
@@ -6,12 +6,7 @@
#include <Sketch.hpp>
template
<
size_t
Dimension
>
class
DiscreteBoundary
{
};
template
<
>
class
DiscreteBoundary
<
2
>
{
public:
DiscreteBoundary
()
{};
...
...
src/Physics/src/DiscreteRegion.h
View file @
01d8f973
...
...
@@ -7,12 +7,8 @@
#include "MaterialProperties.h"
template
<
size_t
Dimension
>
class
DiscreteRegion
{
};
template
<
>
class
DiscreteRegion
<
2
>
{
class
DiscreteRegion
{
public:
DiscreteRegion
()
:
Region
(
nullptr
),
Material
{
Air
()}
{};
...
...
src/Physics/src/FiniteElementMesh.h
View file @
01d8f973
...
...
@@ -16,7 +16,7 @@ class FiniteElementMeshInterface {
public:
FiniteElementMeshInterface
()
{};
FiniteElementMeshInterface
(
std
::
vector
<
XY
>
nodes
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
r
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b
)
:
Nodes
(
nodes
),
Regions
(
r
),
Boundaries
(
b
)
{
FiniteElementMeshInterface
(
std
::
vector
<
XY
>
nodes
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
r
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b
)
:
Nodes
(
nodes
),
Regions
(
r
),
Boundaries
(
b
)
{
sort_boundaries
();
sort_regions
();
};
...
...
@@ -29,25 +29,25 @@ public:
auto
&
boundary
(
size_t
i
)
{
return
Boundaries
[
i
];
};
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
boundary
(
std
::
shared_ptr
<
Curve
const
>
const
&
c
)
const
{
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
boundary
(
std
::
shared_ptr
<
Curve
const
>
const
&
c
)
const
{
// Boundary may be discontinuous (e.g. multiple overlapping curves). Therefore, multiple DiscreteBoundaries may be returned (in general, upper != lower)
auto
lower_comp
=
[](
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
const
&
x
,
size_t
const
&
y
)
{
return
(
size_t
)
(
x
->
curve
().
get
())
<
y
;
};
auto
upper_comp
=
[](
size_t
const
&
y
,
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
const
&
x
)
{
return
y
<
(
size_t
)
(
x
->
curve
().
get
());
};
auto
lower_comp
=
[](
std
::
shared_ptr
<
DiscreteBoundary
>
const
&
x
,
size_t
const
&
y
)
{
return
(
size_t
)
(
x
->
curve
().
get
())
<
y
;
};
auto
upper_comp
=
[](
size_t
const
&
y
,
std
::
shared_ptr
<
DiscreteBoundary
>
const
&
x
)
{
return
y
<
(
size_t
)
(
x
->
curve
().
get
());
};
auto
lower
=
std
::
lower_bound
(
Boundaries
.
begin
(),
Boundaries
.
end
(),
(
size_t
)
c
.
get
(),
lower_comp
);
auto
upper
=
std
::
upper_bound
(
Boundaries
.
begin
(),
Boundaries
.
end
(),
(
size_t
)
c
.
get
(),
upper_comp
);
if
(
lower
!=
Boundaries
.
end
())
{
return
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
(
lower
,
upper
);
return
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
(
lower
,
upper
);
}
else
{
return
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
();
return
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
();
}
};
// TODO: virtual std::shared_ptr<DiscreteRegion
<2>
> region(std::shared_ptr<Region const> const &r) const = 0;
// TODO: virtual std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Region const> const &r) const = 0;
virtual
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
make_discontinuous
(
std
::
shared_ptr
<
Curve
const
>
const
&
c
)
=
0
;
virtual
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
make_discontinuous
(
std
::
shared_ptr
<
Curve
const
>
const
&
c
)
=
0
;
auto
const
&
boundaries
()
const
{
return
Boundaries
;
};
...
...
@@ -63,15 +63,15 @@ public:
auto
const
&
region
(
size_t
i
)
const
{
return
Regions
[
i
];
};
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>
region
(
std
::
shared_ptr
<
Contour
const
>
const
&
c
)
const
{
auto
comp
=
[](
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>
const
&
x
,
size_t
const
&
y
)
{
return
(
size_t
)
(
x
->
region
().
get
())
<
y
;
};
std
::
shared_ptr
<
DiscreteRegion
>
region
(
std
::
shared_ptr
<
Contour
const
>
const
&
c
)
const
{
auto
comp
=
[](
std
::
shared_ptr
<
DiscreteRegion
>
const
&
x
,
size_t
const
&
y
)
{
return
(
size_t
)
(
x
->
region
().
get
())
<
y
;
};
auto
iter
=
std
::
lower_bound
(
Regions
.
begin
(),
Regions
.
end
(),
(
size_t
)
c
.
get
(),
comp
);
if
(
iter
!=
Regions
.
end
())
{
return
*
iter
;
}
else
{
return
std
::
make_shared
<
DiscreteRegion
<
2
>
>
();
return
std
::
make_shared
<
DiscreteRegion
>
();
}
}
...
...
@@ -91,19 +91,29 @@ public:
virtual
std
::
vector
<
std
::
vector
<
XY
>>
quadrature_points
()
const
=
0
;
virtual
size_t
element_order
()
const
=
0
;
virtual
size_t
nodes_per_element
()
const
=
0
;
virtual
size_t
quadrature_order
()
const
=
0
;
virtual
size_t
quadrature_size
()
const
=
0
;;
virtual
double_t
quadrature_weight
(
size_t
i
)
const
=
0
;
protected:
std
::
vector
<
XY
>
Nodes
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
Regions
;
// Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
Boundaries
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
Regions
;
// Contains vector of size_t referencing Triangles (and later Quadrilaterals)
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
Boundaries
;
void
sort_boundaries
()
{
auto
comp
=
[](
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
const
&
x
,
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
const
&
y
)
{
return
(
size_t
)
(
x
->
curve
().
get
())
<
(
size_t
)
(
y
->
curve
().
get
());
};
auto
comp
=
[](
std
::
shared_ptr
<
DiscreteBoundary
>
const
&
x
,
std
::
shared_ptr
<
DiscreteBoundary
>
const
&
y
)
{
return
(
size_t
)
(
x
->
curve
().
get
())
<
(
size_t
)
(
y
->
curve
().
get
());
};
std
::
sort
(
Boundaries
.
begin
(),
Boundaries
.
end
(),
comp
);
}
void
sort_regions
()
{
auto
comp
=
[](
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>
const
&
x
,
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>
const
&
y
)
{
return
(
size_t
)
(
x
->
region
().
get
())
<
(
size_t
)
(
y
->
region
().
get
());
};
auto
comp
=
[](
std
::
shared_ptr
<
DiscreteRegion
>
const
&
x
,
std
::
shared_ptr
<
DiscreteRegion
>
const
&
y
)
{
return
(
size_t
)
(
x
->
region
().
get
())
<
(
size_t
)
(
y
->
region
().
get
());
};
std
::
sort
(
Regions
.
begin
(),
Regions
.
end
(),
comp
);
}
};
...
...
@@ -113,7 +123,7 @@ class FiniteElementMesh : public FiniteElementMeshInterface {
public:
FiniteElementMesh
()
{};
FiniteElementMesh
(
std
::
vector
<
XY
>
nodes
,
std
::
vector
<
Triangle
<
ElementOrder
>>
tris
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
r
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b
)
:
FiniteElementMeshInterface
{
nodes
,
r
,
b
},
Triangles
(
tris
)
{};
FiniteElementMesh
(
std
::
vector
<
XY
>
nodes
,
std
::
vector
<
Triangle
<
ElementOrder
>>
tris
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
r
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b
)
:
FiniteElementMeshInterface
{
nodes
,
r
,
b
},
Triangles
(
tris
)
{};
FiniteElementMesh
(
Mesh
&
m
)
{
// std::vector<XY> nodes
...
...
@@ -121,7 +131,7 @@ public:
// std::vector<std::shared_ptr<Region<2>>> r
// std::vector<std::shared_ptr<Boundary<2>>> b
Nodes
.
reserve
(
m
.
size_points
());
// TODO: This will be larger when Order > 1
Nodes
.
reserve
(
m
.
size_points
());
// TODO: This will be larger when
Element
Order > 1
for
(
size_t
i
=
0
;
i
!=
m
.
size_points
();
++
i
)
{
Nodes
.
emplace_back
(
m
.
point
(
i
));
}
...
...
@@ -144,12 +154,12 @@ public:
Boundaries
.
reserve
(
m
.
size_boundary_constraints
());
for
(
size_t
i
=
0
;
i
!=
m
.
size_boundary_constraints
();
++
i
)
{
Boundaries
.
push_back
(
std
::
make_shared
<
DiscreteBoundary
<
2
>
>
(
m
.
curve
(
i
),
m
.
boundary_nodes
(
i
)));
Boundaries
.
push_back
(
std
::
make_shared
<
DiscreteBoundary
>
(
m
.
curve
(
i
),
m
.
boundary_nodes
(
i
)));
}
Regions
.
reserve
(
m
.
size_contours
());
for
(
size_t
i
=
0
;
i
!=
m
.
size_contours
();
++
i
)
{
Regions
.
push_back
(
std
::
make_shared
<
DiscreteRegion
<
2
>
>
(
m
.
contour
(
i
)));
// TODO: Assign material properties here
Regions
.
push_back
(
std
::
make_shared
<
DiscreteRegion
>
(
m
.
contour
(
i
)));
// TODO: Assign material properties here
}
Triangles
.
reserve
(
m
.
size_triangles
());
...
...
@@ -169,6 +179,16 @@ public:
size_t
elements_size
()
const
override
{
return
Triangles
.
size
();
};
size_t
element_order
()
const
override
{
return
ElementOrder
;
};
size_t
nodes_per_element
()
const
override
{
return
Triangle
<
ElementOrder
>::
NumNodes
;
};
size_t
quadrature_order
()
const
override
{
return
QuadratureOrder
;
};
size_t
quadrature_size
()
const
override
{
return
TriangleQuadrature
<
QuadratureOrder
>::
size
;
};
double_t
quadrature_weight
(
size_t
i
)
const
override
{
return
TriangleQuadrature
<
QuadratureOrder
>::
w
[
i
];
};
DiagonalMatrixGroup
determinant
()
const
override
{
DiagonalMatrixGroup
mat
(
TriangleQuadrature
<
QuadratureOrder
>::
size
,
Triangles
.
size
());
...
...
@@ -210,16 +230,16 @@ public:
return
qp
;
}
// TODO: std::shared_ptr<DiscreteRegion
<2>
> region(std::shared_ptr<Region const> const &r) const override {...}
// TODO: std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Region const> const &r) const override {...}
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
make_discontinuous
(
std
::
shared_ptr
<
Curve
const
>
const
&
c
)
override
{
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
make_discontinuous
(
std
::
shared_ptr
<
Curve
const
>
const
&
c
)
override
{
// TODO: Should this be implemented in the refineable mesh instead? It might be simpler.
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b_vec
=
boundary
(
c
);
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b_vec
=
boundary
(
c
);
if
(
b_vec
.
size
()
==
0
||
b_vec
.
size
()
==
2
)
{
// Boundary curve not found || boundary pair found
return
b_vec
;
}
else
if
(
b_vec
.
size
()
>
2
)
{
// should not happen
throw
std
::
runtime_error
(
"FiniteElementMesh<2>::make_discontinuous expects to find at most 2 DiscreteBoundary
<2>
objects"
);
throw
std
::
runtime_error
(
"FiniteElementMesh<2>::make_discontinuous expects to find at most 2 DiscreteBoundary objects"
);
}
// else b.size() == 1 and we need to construct an new boundary
// Copy nodes
...
...
@@ -233,7 +253,7 @@ public:
// Construct new DiscontinuousBoundary
bool
is_closed
=
b_vec
[
0
]
->
is_closed
();
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
db
=
std
::
make_shared
<
DiscreteBoundary
<
2
>
>
(
b_vec
[
0
]
->
curve
(),
discontinuous_nodes
,
is_closed
);
std
::
shared_ptr
<
DiscreteBoundary
>
db
=
std
::
make_shared
<
DiscreteBoundary
>
(
b_vec
[
0
]
->
curve
(),
discontinuous_nodes
,
is_closed
);
Boundaries
.
push_back
(
db
);
sort_boundaries
();
...
...
@@ -249,7 +269,7 @@ public:
double_t
dy0
=
Nodes
[
boundary_nodes
[
ii
]].
y
()
-
Nodes
[
boundary_nodes
[
i
]].
y
();
for
(
Triangle
<
ElementOrder
>
&
t
:
Triangles
)
{
for
(
size_t
k
=
0
;
k
!=
3
;
++
k
)
{
for
(
size_t
k
=
0
;
k
!=
3
;
++
k
)
{
// TODO: This essentially assumes Element order = 1
if
(
t
.
node
(
k
)
==
boundary_nodes
[
i
])
{
size_t
k0
=
t
.
node
(
0
);
size_t
k1
=
t
.
node
(
1
);
...
...
src/Physics/src/Forcing.h
View file @
01d8f973
...
...
@@ -19,15 +19,10 @@ public:
virtual
Eigen
::
VectorXd
operator
()(
double
const
t
)
const
=
0
;
};
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
HomogeneousForcing
:
public
Forcing
{
};
template
<
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
HomogeneousForcing
<
2
,
ElementOrder
,
QuadratureOrder
>
:
public
Forcing
{
public:
HomogeneousForcing
(
std
::
function
<
double
(
double
)
>
f
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
regions
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
regions
,
SparseMatrixGroup
const
&
basis
,
DiagonalMatrixGroup
const
&
weights
)
:
Function
{
f
},
Regions
{
regions
},
Basis
(
basis
),
Weights
(
weights
)
{
...
...
@@ -44,7 +39,7 @@ public:
Eigen
::
VectorXd
operator
()(
double
const
t
)
const
override
{
Eigen
::
VectorXd
output
=
(
Basis
[
0
]
*
(
Weights
(
0
).
matrix
().
asDiagonal
()
*
Indicator
));
// TODO: ?SparseVector?
for
(
size_t
i
=
1
;
i
!=
TriangleQuadrature
<
QuadratureOrder
>::
size
;
++
i
)
{
for
(
size_t
i
=
1
;
i
!=
Basis
.
size
()
;
++
i
)
{
output
+=
(
Basis
[
i
]
*
(
Weights
(
i
).
matrix
().
asDiagonal
()
*
Indicator
));
}
...
...
@@ -57,7 +52,7 @@ public:
protected:
std
::
function
<
double
(
double
)
>
Function
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
Regions
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
Regions
;
SparseMatrixGroup
const
&
Basis
;
DiagonalMatrixGroup
const
&
Weights
;
...
...
src/Physics/src/MatrixGroup.h
View file @
01d8f973
...
...
@@ -4,6 +4,8 @@
#include "Eigen"
#include "Eigen/Sparse"
// TODO: MatrixGroup interface class
class
SparseMatrixGroup
{
public:
SparseMatrixGroup
(
size_t
const
Q
,
size_t
const
rows
,
size_t
const
cols
,
size_t
const
nnz
)
{
...
...
@@ -20,6 +22,8 @@ public:
auto
&
operator
[](
size_t
const
q
)
const
{
return
Matrices
[
q
];
};
size_t
size
()
const
{
return
Matrices
.
size
();
};
void
transform
(
Eigen
::
SparseMatrix
<
double
>
&
A
)
{
for
(
size_t
i
=
0
;
i
!=
Matrices
.
size
();
++
i
)
{
Matrices
[
i
]
=
A
*
Matrices
[
i
];
...
...
@@ -40,8 +44,11 @@ public:
}
auto
&
operator
()(
size_t
const
q
)
const
{
return
Matrices
[
q
];
};
auto
&
operator
()(
size_t
const
q
)
{
return
Matrices
[
q
];
};
size_t
size
()
const
{
return
Matrices
.
size
();
};
protected:
std
::
vector
<
Eigen
::
ArrayXd
>
Matrices
;
};
...
...
@@ -58,6 +65,8 @@ public:
auto
&
dy
(
size_t
const
q
,
size_t
const
i
,
size_t
const
j
)
{
return
Dy
(
q
,
i
,
j
);
};
size_t
size
()
const
{
return
Dx
.
size
();
};
void
transform
(
Eigen
::
SparseMatrix
<
double
>
&
A
)
{
Dx
.
transform
(
A
);
Dy
.
transform
(
A
);
...
...
src/Physics/src/Physics.h
View file @
01d8f973
...
...
@@ -19,31 +19,9 @@
// TODO: Higher level classes should implement interfaces with appropriate physics grammar
// TODO: Implement material property dependencies using pointers to virtual member functions
enum
class
FieldVariable
{
MagneticVectorPotential
=
1
,
A
=
1
,
MagneticFluxDensity
=
2
,
B
=
2
,
ElectricFluxDensity
=
4
,
D
=
4
,
ElectricFieldIntensity
=
5
,
E
=
5
,
MagneticFieldIntensity
=
8
,
H
=
8
,
ElectricCurrentDensity
=
10
,
J
=
10
,
ElectricScalarPotential
=
26
+
21
,
Phi
=
26
+
21
,
MagneticScalarPotential
=
26
+
23
,
Psi
=
26
+
23
};
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
PhysicsInterface
{
public:
static
constexpr
size_t
QuadraturePoints
=
TriangleQuadrature
<
QuadratureOrder
>::
size
;
/* TODO: template<size_t Dimension, size_t ElementOrder, size_t QuadratureOrder> class PhysicsData : public PhysicsInterface {...};
/*
virtual Eigen::VectorXd init_residual() = 0;
virtual Eigen::SparseMatrix<double> init_jacobian() = 0;
virtual void residual(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &r) = 0;
...
...
@@ -51,11 +29,14 @@ public:
virtual void apply_conditions() = 0;
*/
//PhysicsInterface(FiniteElementMesh<ElementOrder, QuadratureOrder> &d) : Domain{d},
PhysicsInterface
(
std
::
shared_ptr
<
FiniteElementMeshInterface
>
d
)
:
Domain
{
d
},
ElementWeights
{
QuadraturePoints
,
d
->
elements_size
()},
Derivatives
{
QuadraturePoints
,
d
->
nodes_size
(),
d
->
elements_size
(),
Triangle
<
ElementOrder
>::
NumNodes
},
Basis
{
QuadraturePoints
,
d
->
nodes_size
(),
d
->
elements_size
(),
Triangle
<
ElementOrder
>::
NumNodes
}
{}
ElementWeights
{
d
->
quadrature_size
(),
d
->
elements_size
()},
Derivatives
{
d
->
quadrature_size
(),
d
->
nodes_size
(),
d
->
elements_size
(),
d
->
nodes_per_element
()},
Basis
{
d
->
quadrature_size
(),
d
->
nodes_size
(),
d
->
elements_size
(),
d
->
nodes_per_element
()}
{
std
::
cout
<<
"TODO: Think about storing matrices in the FiniteElementMesh object since this class is simply forwarding it's properties"
<<
std
::
endl
;
std
::
cout
<<
"TODO: Rather, give FiniteElementMesh methods which produce ElementWeight, Derivative, and Basis matrices of appropriate sizes"
<<
std
::
endl
;
std
::
cout
<<
"TODO: That is, FiniteElementMesh is essentially a MatrixFactory"
<<
std
::
endl
;
};
double
time
()
const
{
return
Time
;
};
...
...
@@ -63,8 +44,6 @@ public:
void
time
(
double
t
)
{
Time
=
t
;
};
//FiniteElementMesh<ElementOrder, QuadratureOrder> const &domain() const { return Domain; };
std
::
shared_ptr
<
FiniteElementMeshInterface
>
domain
()
const
{
return
Domain
;
};
DiagonalMatrixGroup
const
&
weights
()
{
return
ElementWeights
;
};
...
...
@@ -78,8 +57,6 @@ public:
protected:
double
Time
{
0.0
};
//FiniteElementMesh<ElementOrder, QuadratureOrder> &Domain;
std
::
shared_ptr
<
FiniteElementMeshInterface
>
Domain
;
std
::
vector
<
std
::
shared_ptr
<
Forcing
>>
ForcingCondtions
;
...
...
@@ -92,16 +69,9 @@ protected:
Eigen
::
SparseMatrix
<
double
>
BoundaryConditionMatrix
;
};
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
,
FieldVariable
V
>
class
Magnetostatic
:
public
PhysicsInterface
<
Dimension
,
ElementOrder
,
QuadratureOrder
>
{
};
template
<
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
Magnetostatic
<
2
,
ElementOrder
,
QuadratureOrder
,
FieldVariable
::
MagneticVectorPotential
>
:
public
PhysicsInterface
<
2
,
ElementOrder
,
QuadratureOrder
>
{
class
Magnetostatic
:
public
PhysicsInterface
{
public:
//Magnetostatic(FiniteElementMesh<ElementOrder, QuadratureOrder> &d) : PhysicsInterface<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.size_nodes()}, Elements{d.elements_size()} {}
Magnetostatic
(
std
::
shared_ptr
<
FiniteElementMeshInterface
>
d
)
:
PhysicsInterface
<
2
,
ElementOrder
,
QuadratureOrder
>
{
d
},
Unknowns
{
d
->
nodes_size
()},
Elements
{
d
->
elements_size
()}
{}
Magnetostatic
(
std
::
shared_ptr
<
FiniteElementMeshInterface
>
d
)
:
PhysicsInterface
{
d
},
Unknowns
{
d
->
nodes_size
()},
Elements
{
d
->
elements_size
()}
{}
Eigen
::
SparseMatrix
<
double
>
init_unknown_matrix
()
const
{
return
Eigen
::
SparseMatrix
<
double
>
(
Unknowns
,
Unknowns
);
}
...
...
@@ -179,7 +149,7 @@ public:
r
=
-
f
;
J
.
setZero
();
for
(
size_t
i
=
0
;
i
!=
TriangleQuadrature
<
Q
uadrature
Order
>::
size
;
++
i
)
{
for
(
size_t
i
=
0
;
i
!=
Domain
->
q
uadrature
_
size
()
;
++
i
)
{
// Caclulate flux density
Fx
=
this
->
Derivatives
.
dy
(
i
).
transpose
()
*
v
;
Fy
=
-
this
->
Derivatives
.
dx
(
i
).
transpose
()
*
v
;
...
...
@@ -207,8 +177,8 @@ public:
this
->
Basis
=
this
->
Domain
->
basis
();
this
->
Derivatives
=
this
->
Domain
->
derivative
();
for
(
size_t
i
=
0
;
i
!=
TriangleQuadrature
<
Q
uadrature
Order
>::
size
;
++
i
)
{
this
->
ElementWeights
(
i
)
*=
TriangleQuadrature
<
Q
uadrature
O
rder
>::
w
[
i
];
for
(
size_t
i
=
0
;
i
!=
Domain
->
q
uadrature
_
size
()
;
++
i
)
{
this
->
ElementWeights
(
i
)
*=
Domain
->
quadrature_weight
(
i
);
//
TriangleQuadrature<
(Domain->q
uadrature
_o
rder
())
>::w[i];
}
}
...
...
@@ -220,60 +190,60 @@ public:
}
void
add_current_density
(
std
::
function
<
double
(
double
)
>
func
,
std
::
vector
<
std
::
shared_ptr
<
Contour
const
>>
contours
)
{
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
regions
;
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
regions
;
for
(
auto
&
c
:
contours
)
{
regions
.
push_back
(
this
->
Domain
->
region
(
c
));
}
this
->
ForcingCondtions
.
push_back
(
std
::
make_shared
<
HomogeneousForcing
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
func
,
regions
,
this
->
Basis
,
this
->
ElementWeights
));
this
->
ForcingCondtions
.
push_back
(
std
::
make_shared
<
HomogeneousForcing
>
(
func
,
regions
,
this
->
Basis
,
this
->
ElementWeights
));
}
void
add_current_density
(
std
::
function
<
double
(
double
)
>
func
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
<
2
>
>>
regions
)
{
this
->
ForcingCondtions
.
push_back
(
std
::
make_shared
<
HomogeneousForcing
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
func
,
regions
,
this
->
Basis
,
this
->
ElementWeights
));
void
add_current_density
(
std
::
function
<
double
(
double
)
>
func
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteRegion
>>
regions
)
{
this
->
ForcingCondtions
.
push_back
(
std
::
make_shared
<
HomogeneousForcing
>
(
func
,
regions
,
this
->
Basis
,
this
->
ElementWeights
));
}
void
remove_current_density
(
size_t
i
)
{
this
->
ForcingCondtions
.
erase
(
this
->
ForcingCondtions
.
begin
()
+
i
);
}
void
add_magnetic_insulation
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
boundaries
)
{
// TODO: should boundaries be unique_ptrs?
this
->
BoundaryConditions
.
push_back
(
std
::
make_shared
<
ZeroDirichlet
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
boundaries
));
void
add_magnetic_insulation
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
boundaries
)
{
// TODO: should boundaries be unique_ptrs?
this
->
BoundaryConditions
.
push_back
(
std
::
make_shared
<
ZeroDirichlet
>
(
boundaries
));
}
void
add_magnetic_insulation
(
std
::
vector
<
std
::
shared_ptr
<
Curve
const
>>
boundaries
)
{
this
->
BoundaryConditions
.
push_back
(
std
::
make_shared
<
ZeroDirichlet
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
boundaries
,
this
->
Domain
));
this
->
BoundaryConditions
.
push_back
(
std
::
make_shared
<
ZeroDirichlet
>
(
boundaries
,
this
->
Domain
));
}
// TODO: Using forwarding to support many constructors
auto
add_periodic_boundary
(
std
::
vector
<
VariableMap
>
map
,
bool
antiperiodic
)
{
auto
pbc
=
std
::
make_shared
<
PeriodicBoundaryCondition
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
map
,
antiperiodic
);
auto
pbc
=
std
::
make_shared
<
PeriodicBoundaryCondition
>
(
map
,
antiperiodic
);
this
->
BoundaryConditions
.
push_back
(
pbc
);
return
pbc
;
}
auto
add_periodic_boundary
(
std
::
vector
<
PeriodicBoundaryPair
>
periodic_boundary_pairs
,
bool
antiperiodic
)
{
auto
pbc
=
std
::
make_shared
<
PeriodicBoundaryCondition
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
periodic_boundary_pairs
,
this
->
Domain
,
antiperiodic
);
auto
pbc
=
std
::
make_shared
<
PeriodicBoundaryCondition
>
(
periodic_boundary_pairs
,
this
->
Domain
,
antiperiodic
);
this
->
BoundaryConditions
.
push_back
(
pbc
);
return
pbc
;
}
auto
add_periodic_boundary
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b0
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>>
b1
,
std
::
vector
<
bool
>
orientation
,
auto
add_periodic_boundary
(
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b0
,
std
::
vector
<
std
::
shared_ptr
<
DiscreteBoundary
>>
b1
,
std
::
vector
<
bool
>
orientation
,
bool
antiperiodic
)
{
auto
pbc
=
std
::
make_shared
<
PeriodicBoundaryCondition
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
b0
,
b1
,
orientation
,
antiperiodic
);
auto
pbc
=
std
::
make_shared
<
PeriodicBoundaryCondition
>
(
b0
,
b1
,
orientation
,
antiperiodic
);
this
->
BoundaryConditions
.
push_back
(
pbc
);
return
pbc
;
}
auto
add_sliding_interface
(
std
::
shared_ptr
<
Curve
const
>
interface
,
bool
antiperiodic
)
{
auto
boundary
=
this
->
Domain
->
make_discontinuous
(
interface
);
auto
condition
=
std
::
make_shared
<
SlidingInterface
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
boundary
[
0
],
boundary
[
1
],
antiperiodic
);
auto
condition
=
std
::
make_shared
<
SlidingInterface
>
(
boundary
[
0
],
boundary
[
1
],
antiperiodic
);
this
->
BoundaryConditions
.
push_back
(
condition
);
return
condition
;
}
auto
add_sliding_interface
(
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
b0
,
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
b1
,
bool
antiperiodic
)
{
auto
condition
=
std
::
make_shared
<
SlidingInterface
<
2
,
ElementOrder
,
QuadratureOrder
>
>
(
b0
,
b1
,
antiperiodic
);
auto
add_sliding_interface
(
std
::
shared_ptr
<
DiscreteBoundary
>
b0
,
std
::
shared_ptr
<
DiscreteBoundary
>
b1
,
bool
antiperiodic
)
{
auto
condition
=
std
::
make_shared
<
SlidingInterface
>
(
b0
,
b1
,
antiperiodic
);
this
->
BoundaryConditions
.
push_back
(
condition
);
return
condition
;
}
...
...
src/Physics/src/PostProcessor.h
View file @
01d8f973
...
...
@@ -3,12 +3,11 @@
#include "Physics.h"
template
<
size_t
Dimension
,
size_t
ElementOrder
,
size_t
QuadratureOrder
>
class
PostProcessor
{
public:
PostProcessor
(
std
::
shared_ptr
<
PhysicsInterface
<
Dimension
,
ElementOrder
,
QuadratureOrder
>
>
p
)
:
Physics
{
p
}
{}
PostProcessor
(
std
::
shared_ptr
<
PhysicsInterface
>
p
)
:
Physics
{
p
}
{}
std
::
shared_ptr
<
PhysicsInterface
<
Dimension
,
ElementOrder
,
QuadratureOrder
>
>
Physics
;
std
::
shared_ptr
<
PhysicsInterface
>
Physics
;
};
#endif //OERSTED_POSTPROCESSOR_H
test/LibraryIntegration/Mesh_To_FEM_Test.cpp
View file @
01d8f973
...
...
@@ -40,14 +40,14 @@ TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
me
.
refine
();
me
.
save_as
(
SAVE_DIR
,
std
::
string
(
"ful
L
_circle_mesh"
));
me
.
save_as
(
SAVE_DIR
,
std
::
string
(
"ful
l
_circle_mesh"
));
// Convert to FiniteElementMesh
//FiniteElementMesh<1, 1> fem{me};
std
::
shared_ptr
<
FiniteElementMeshInterface
>
fem
=
std
::
make_shared
<
FiniteElementMesh
<
1
,
1
>>
(
me
);
for
(
std
::
shared_ptr
<
DiscreteBoundary
<
2
>
>
const
&
b
:
fem
->
boundaries
())
{
for
(
std
::
shared_ptr
<
DiscreteBoundary
>
const
&
b
:
fem
->
boundaries
())
{