Filetype PDF | Posted on 24 Jan 2023 | 2 years ago
Authentication
digital resources
117x Filetype PDF File size 0.21 MB Source: dune-project.org
{ 2, 3 } }; // S :=2,S :=3 gt.makeSimplex(2); // same as makeTriangle() (0,0,1)
10 11 3 (0,0,1)
In General FieldMatrix Q(k); // Qij :=k ∀i,j gt.makeCube(3); // same as makeHexahedron() 3
// for each makeShape() there is an isShape()
main() template assert(i < S.rows()); // get number of rows assert(gt.isHexahedron());
assert(j < S.cols()); // get number of columns assert(gt.isCube()); // ignore dimension 5
#include k = S[i][j]; // access entry assert(gt.dim() == 3); // check dimension
S[i][j] = k; // assign entry 3 3 4
for(const auto &row : S) Concept Geometry right
#include 2 2 2
for(const auto &entry : row) left (0,1,0)
using Geo = ...; Geo geo; 1 (0,1,0)
int main(int argc, char **argv) k += entry; // access each entry front 1 2
for(auto &row : S) 0
{ using ctype = Geo::ctype; bottom
Dune::MPIHelper::instance(argc, argv); for(auto &entry : row) 0 1 0 0 1
entry = k; // modify each entry int ldim = Geo::mydimension; // local dim (0,0,0) (1,0,0)(0,0,0) (1,0,0)
int gdim = Geo::coorddimension; // global dim (0,0,1)
(0,0,1)
// your code goes here auto L = Q.leftmultiplyany(S); // L:=SQ 4 4
} auto R = Q.rightmultiplyany(S); // R:=QS Geo::LocalCoordinate xl; // xˆ ∈ ctypeldim
Q.leftmultiply(S); // S := SQ Geo::GlobalCoordinate x; // x∈ctypegdim
.cc-file template x = geo.global(xl); // x:=g(xˆ) 4 6 7
Q.rightmultiply(S); // S := QS −1 rear
#include S += Q; S -= Q; // S :=S+Q, S :=S−Q xl = geo.local(x); // xˆ := g (x) 1 2 3 4 5
left 2 2 3 3
(0,1,0) (1,1,0)
S.axpy(k, Q); // S := S +kQ p 3 (0,1,0) (1,1,0)
−1 // J−T ∈ctypegdim×ldim, J := ∂g /∂xˆ , µ := | detJTJ| front 0 right 0 1
// your code and includes go here S *= k; S /= k; // S :=kS, S :=k S ij i j bottom
S.invert(); // S := S−1 Geo::JacobianInverseTransposed JInvT = 0 1 0 2 1
.hh-file template q geo.jacobianInverseTransposed(xl); (0,0,0) (1,0,0) (0,0,0) (1,0,0)
P (0,1,1)
2 ctype mu = geo.integrationElement(xl); (0,1,1)
r = S.frobenius_norm(); // r:= ij |Sij|
// For a header that is included like r = S.infinity_norm(); // r := max {P |S |} 5
// #include i j ij 5
k = S.determinant(); // k := detS GeometryType gt = geo.type(); // shape 4
#ifndef DUNE_MODULE_DIR_HEADER_NAME_HH assert(i < geo.corners()); // count corners (0,0,1) (1,0,1)
#define DUNE_MODULE_DIR_HEADER_NAME_HH top 7 8
Q.mv (x, y); // y :=Qx x = geo.corner(i); // access corner 3 4 (0,0,1) (1,0,1)
Q.mtv (x, y); // y :=QTx x = geo.center(); // roughly 3 6 4
// your code and includes go here Q.umv (x, y); // y :=y+Qx ctype v = geo.volume(); // in global coords 2
// do not #include Q.umtv (x, y); // y :=y+QTx 2
H template 1 right
#endif // DUNE_MODULE_DIR_HEADER_NAME_HH Q.umhv (x, y); // y :=y+Q x class ReferenceElements;
Q.usmv (k, x, y); // y :=y+kQx left 0
T front 0 1
Q.usmtv(k, x, y); // y :=y+kQ x #include 2 2
Q.usmhv(k, x, y); // y :=y+kQHx (0,1,0)
❞✉♥❡✲❝♦♠♠♦♥ (0,1,0)
Q.solve (x, y); // find x such that Qx=y using Factory = ReferenceElements; 3 4 5
In the following, r is of type R, which may be a scalar real type, const auto &refTet = Factory::simplex(); bottom
#define DUNE_THROW(ExceptionType, message) const auto &refHex = Factory::cube(); 0 1 0 3 1
e.g. double or float. k is of type K, which may be may be R or GeometryType gt; gt.makePrism(); (0,0,0) (1,0,0)(0,0,0) (1,0,0)
std::complex. #include (0,1,1) (1,1,1)
const auto &ref = Factory::general(gt); 6 7 (0,1,1) (1,1,1)
template class FieldVector; 6 11 7
if(i > limit) 5
// Info about ref itself (0,0,1) top (1,0,1) 8 9
#include DUNE_THROW(Exception, "Error: i > limit (" (0,0,1) (1,0,1)
gt = ref.type(); 4 5 4 10 5
<< i << " > " << limit << ")"); ctype v = ref.volume(); 3
FieldVector x = { 0, 1 }; // x :=0,x :=1 rear 2 3
0 1 ref.size(c); // count subentities of codim c
FieldVector y(k); // x :=k ∀i template std::string className();
i 0 1
template std::string className(T& t); left 2 right 0 1
assert(i < x.dim()); // get number of entries // Info about subentity (i,c) front
#include gt = ref.type(i,c); 2 3 2 7 3
k = x[i]; x[i] = k; // access/assign entry // position of barycenter (0,1,0) (1,1,0) (0,1,0) (1,1,0)
for(const auto &entry : x) 4 4 5
template FieldVector x = ref.position(i,c) bottom
k += entry; // access each entry // count sub-subentities of codim cc 0 1 0 6 1
for(auto &entry : x) void printTypes(const Vector &v) { (0,0,0) (1,0,0) (0,0,0) (1,0,0)
entry = k; // modify each entry std::cerr << "Info: Vector type is " ref.size(i,c, cc);
<< className() // transform number of sub-subentity to ref template class QuadratureRules;
x += y; x -= y; // x:=x+y, x:=x−y << ", entry type is " ref.subEntity(i,c, ii,cc);
x *= k; x /= k; // x:=kx, x:=k−1x << className(v[0]) << std::endl; (0,1) (0,1) (1,1) #include
k = x * y; // k := xTy = x·y = P x y } 3
i i i 2 2 3
k = x.dot(y); // k := xHy = x¯·y = P x¯ y K f(const FieldVector &x);
i i i int p; // max polynomial order of f
r = x.one_norm(); // r := P |xi| ❞✉♥❡✲❣❡♦♠❡tr②
pi
r = x.two_norm(); // r := P|xi|2 K result = 0;
i class GeometryType; GeometryType gt; gt.makeSimplex(dim);
r = x.infinity_norm(); // r:=maxi{|xi|} 1 2 0 1 for(const auto &qp :
template #include QuadratureRules::rule(gt, p))
class FieldMatrix; result += qp.weight() * f(qp.position());
GeometryType gt; // now result contains the integral of f()
#include gt.makeVertex(); gt.makeLine(); // over the reference-simplex of dimension dim
gt.makeTriangle(); gt.makeQuadrilateral(); 0 0 1 0 2 1 TODO:integral over a geometry over a scalar
FieldMatrix S = gt.makeTetrahdron(); gt.makePyramid(); (0,0) (1,0) (0,0) (1,0) TODO: integral over a geometry over a gradient (incl piola)
{ { 0, 1 }, // S :=0,S :=1 gt.makePrism(); gt.makeHexahedron();
00 01
1
int codim = Entity::codimension; class LocalKey; – DoF position info format
int dim = Entity::dimension; // as on Grid // interpolate f to piecewise constants
❞✉♥❡✲❣r✐❞ TODO
int mydim = Entity::mydimension; std::vector p0(gv.size(0)) TODO: list of local finite elements
Grid (YaspGrid, UGGrid, OneDGrid, GeometryGrid) for(const auto &e : elements(gv))
GridView (LevelGridView, LeafGridView) GeometryType gt = e.type(); // Shape p0[set.index(e)] = f(e.geometry().center()); ❞✉♥❡✲✐st❧
IndexSet // the LevelGridView that e is part of
Entity (elements, facets, edges, vertices) int l = e.level(); // interpolate f to P1/Q1
std::vector p1(gv.size(dim)); BlockVector BCRSMatrix
Geometry (entity to global) // transform mydimension -> dimensionworld for(const auto &v : vertices(gv)) MatrixAdapter Preconditioner
Intersection Entity::Geometry geo = e.geometry(); p1[set.index(v)] = f(v.geometry().center()); list of preconditioners
Geometry (intersection to global) Solver Interface InverseOperatorResult
Entity (inside/outside element/cell) Concept Intersection – connectivity between elements // output the two interpolations of f list of solvers
Geometry (intersection to inside/outside) VTKWriter writer(gv);
Concept Grid – hierarchy of meshes Intersection isect; writer.addCellData(p0, "constants");
writer.addVertexData(p1, "linears");
Grid g; using ctype = Intersection::ctype; writer.write("file_name_base");
// local coords (== Grid::dimension - 1)
using ctype = Grid::ctype; int mydim = Intersection::mydimension;
int dim = Grid::dimension; // global coords (== Grid::dimensionworld) ❞✉♥❡✲❧♦❝❛❧❢✉♥❝t✐♦♥s
// think "surface grid" int dimw = Intersection::dimensionworld;
int dimw = Grid::dimensionworld; LocalFiniteElement
GeometryType gt = isect.type(); // Shape
g.globalRefine(n); // add n levels LocalBasis
assert(g.maxLevel() > 0); // transform intersection -> world LocalInterpolation
// all coarse/macro entities Intersection::Geometry geo = isect.geometry(); LocalCoefficients
auto levelView = g.levelGridView(0); Intersection::LocalCoordinate xl; Concept LocalFiniteElement
// all finest/leaf entities Intersection::GlobalCoordinate nu_u, nu_q;
auto leafView = g.leafGridView(); // kν k =1, ν :=ν ·geo.integrationElement(xl) using LocalFiniteElement = ...;
u 2 q u
nu_u = isect.unitOuterNormal(xl); LocalFiniteElement lfe;
Concept GridView – one mesh from the hierarchy nu_q = isect.integrationOuterNormal(xl);
GeometryType gt = lfe.type();
GridView gv; using Element = Intersection::Entity; unsigned shapeFunctionCount = lfe.size();
using LGeo = Intersection::LocalGeometry;
using Grid = GridView::Grid; const auto &lb = lfe.localBasis();
using ctype = GridView::ctype; // as on Grid // inside element (always exists) const auto &li = lfe.localInterpolation();
int dim = GridView::dimension; Element in = isect.inside(); const auto &lc = lfe.localCoefficients();
int dimw = GridView::dimensionworld; // transform intersection -> inside
LGeo inGeo = isect.geometryInInside(); Concept LocalBasis – evaluate shape functions and derivatives
const Grid &g = gv.grid(); // index of subfacet of in that contains isect
const auto &idxSet = gv.indexSet(); int inIdx = isect.indexInInside(); using LocalBasis = ...; LocalBasis lb;
// count entities...
int n = gv.size(c); // with codim c if(isect.neighbor()) { // check outside exists unsigned shapeFunctionCount = lb.size();
int n = gv.size(gt); // with GeometryType gt Element out = isect.outside(); unsigned p = lb.order(); // max polynom order
LGeo outGeo = isect.geometryInOutside(); using T = LocalBasis::Traits;
// iterate over entities in gv int outIdx = isect.indexInOutside();
for(const auto &elem : elements(gv)) ...; } // otherwise on domain boundary unsigned dimD = T::dimDomain;
for(const auto &facet : facets (gv)) ...; using DF = T::DomainFieldType;
dimD
for(const auto &edge : edges (gv)) ...; template class YaspGrid; using Domain = T::DomainType; // DF
for(const auto &vertex : vertices(gv)) ...; Yet Another Structured Parallel Grid unsigned dimR = T::dimRange;
// iterate intersections of elem in gv Implements concept Grid. using RF = T::RangeFieldType;
dimR
for(const auto &isect : using Range = T::RangeType; // RF
intersections(gv, elem)) ...; #include
Domain xl; // xˆ
Concept IndexSet – numbering within GridViews 2 (i)
// construct unit square [0,1] with one element std::vector y; // y[i][j]=y
YaspGrid<2> grid0({ 1, 1 }, { 1, 1 }); j
Entities of different shape (GeometryType) are numbered sepa- (i) (i)
// resizes y to fit; y := ϕˆ (xˆ) ∀i
rately. See MultipleCodimMultipleGeomTypeMapper. j j
// construct cube [−1,1]3 with 8=23 elements lb.evaluateFunction(xl, y);
const IndexSet &idxSet; YaspGrid<3> grid1({ -1, -1, -1 }, { 1, 1, 1 },
{ 2, 2, 2 }); // only guaranteed for T::diffOrder > 0:
Entity e; // any codim using Jacobian = T::JacobianType; // RFdimR×dimD
int i, c; // number/codim of subentity template class VTKWriter; std::vector J; // J[i][j][k]=J(i)
Generate files for paraview jk
idxSet.index(e); // index of e in gv (i) (i)
// resizes J to fit; J := (∂ϕˆ /∂xˆ )| ∀i,j,k
idxSet.subIndex(e, i,c); // index of subentity jk j k xˆ
#include lb.evaluateJacobian(xl, J);
ˆ (i)
Concept Entity – elements, facets, edges, vertices GridView gv; // for scalar bases (dimR==1): J[i][0]=∇ϕˆ
Entity e; double f(FieldVector xg); Concept LocalInterpolation – interpolate into shape functions
TODO
// all entities: mydim + codim == dim // for multiple possible GeometryTypes use Concept LocalCoefficients – DoF position database
// elements: codim == 0; facets: codim == 1 // MultipleCodimMultipleGeomTypeMapper instead TODO
// edges: mydim == 1; vertices: mydim == 0 const auto &set = gv.indexSet();
2
no reviews yet
Please Login to review.
