Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82000602
MeshTri3.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 9, 01:26
Size
18 KB
Mime Type
text/x-c
Expires
Wed, Sep 11, 01:26 (2 d)
Engine
blob
Format
Raw Data
Handle
20643015
Attached To
rGOOSEFEM GooseFEM
MeshTri3.hpp
View Options
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef GOOSEFEM_MESHTRI3_CPP
#define GOOSEFEM_MESHTRI3_CPP
// -------------------------------------------------------------------------------------------------
#include "MeshTri3.h"
// ===================================== GooseFEM::Mesh::Tri3 ======================================
namespace
GooseFEM
{
namespace
Mesh
{
namespace
Tri3
{
// ------------------------------------------ constructor ------------------------------------------
inline
Regular
::
Regular
(
size_t
nelx
,
size_t
nely
,
double
h
)
:
m_h
(
h
),
m_nelx
(
nelx
),
m_nely
(
nely
)
{
assert
(
m_nelx
>=
1
);
assert
(
m_nely
>=
1
);
m_nnode
=
(
m_nelx
+
1
)
*
(
m_nely
+
1
);
m_nelem
=
m_nelx
*
m_nely
*
2
;
}
// -------------------------------------- number of elements ---------------------------------------
inline
size_t
Regular
::
nelem
()
const
{
return
m_nelem
;
}
// ---------------------------------------- number of nodes ----------------------------------------
inline
size_t
Regular
::
nnode
()
const
{
return
m_nnode
;
}
// ---------------------------------- number of nodes per element ----------------------------------
inline
size_t
Regular
::
nne
()
const
{
return
m_nne
;
}
// ------------------------------------- number of dimensions --------------------------------------
inline
size_t
Regular
::
ndim
()
const
{
return
m_ndim
;
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline
xt
::
xtensor
<
double
,
2
>
Regular
::
coor
()
const
{
xt
::
xtensor
<
double
,
2
>
out
=
xt
::
empty
<
double
>
({
m_nnode
,
m_ndim
});
xt
::
xtensor
<
double
,
1
>
x
=
xt
::
linspace
<
double
>
(
0.0
,
m_h
*
static_cast
<
double
>
(
m_nelx
),
m_nelx
+
1
);
xt
::
xtensor
<
double
,
1
>
y
=
xt
::
linspace
<
double
>
(
0.0
,
m_h
*
static_cast
<
double
>
(
m_nely
),
m_nely
+
1
);
size_t
inode
=
0
;
for
(
size_t
iy
=
0
;
iy
<
m_nely
+
1
;
++
iy
)
{
for
(
size_t
ix
=
0
;
ix
<
m_nelx
+
1
;
++
ix
)
{
out
(
inode
,
0
)
=
x
(
ix
);
out
(
inode
,
1
)
=
y
(
iy
);
++
inode
;
}
}
return
out
;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline
xt
::
xtensor
<
size_t
,
2
>
Regular
::
conn
()
const
{
xt
::
xtensor
<
size_t
,
2
>
out
=
xt
::
empty
<
size_t
>
({
m_nelem
,
m_nne
});
size_t
ielem
=
0
;
for
(
size_t
iy
=
0
;
iy
<
m_nely
;
++
iy
)
{
for
(
size_t
ix
=
0
;
ix
<
m_nelx
;
++
ix
)
{
out
(
ielem
,
0
)
=
(
iy
)
*
(
m_nelx
+
1
)
+
(
ix
);
out
(
ielem
,
1
)
=
(
iy
)
*
(
m_nelx
+
1
)
+
(
ix
+
1
);
out
(
ielem
,
2
)
=
(
iy
+
1
)
*
(
m_nelx
+
1
)
+
(
ix
);
++
ielem
;
out
(
ielem
,
0
)
=
(
iy
)
*
(
m_nelx
+
1
)
+
(
ix
+
1
);
out
(
ielem
,
1
)
=
(
iy
+
1
)
*
(
m_nelx
+
1
)
+
(
ix
+
1
);
out
(
ielem
,
2
)
=
(
iy
+
1
)
*
(
m_nelx
+
1
)
+
(
ix
);
++
ielem
;
}
}
return
out
;
}
// ------------------------------ node-numbers along the bottom edge -------------------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesBottomEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nelx
+
1
});
for
(
size_t
ix
=
0
;
ix
<
m_nelx
+
1
;
++
ix
)
out
(
ix
)
=
ix
;
return
out
;
}
// -------------------------------- node-numbers along the top edge --------------------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesTopEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nelx
+
1
});
for
(
size_t
ix
=
0
;
ix
<
m_nelx
+
1
;
++
ix
)
out
(
ix
)
=
ix
+
m_nely
*
(
m_nelx
+
1
);
return
out
;
}
// ------------------------------- node-numbers along the left edge --------------------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesLeftEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nely
+
1
});
for
(
size_t
iy
=
0
;
iy
<
m_nely
+
1
;
++
iy
)
out
(
iy
)
=
iy
*
(
m_nelx
+
1
);
return
out
;
}
// ------------------------------- node-numbers along the right edge -------------------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesRightEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nely
+
1
});
for
(
size_t
iy
=
0
;
iy
<
m_nely
+
1
;
++
iy
)
out
(
iy
)
=
iy
*
(
m_nelx
+
1
)
+
m_nelx
;
return
out
;
}
// ---------------------- node-numbers along the bottom edge, without corners ----------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesBottomOpenEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nelx
-
1
});
for
(
size_t
ix
=
1
;
ix
<
m_nelx
;
++
ix
)
out
(
ix
-
1
)
=
ix
;
return
out
;
}
// ----------------------- node-numbers along the top edge, without corners ------------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesTopOpenEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nelx
-
1
});
for
(
size_t
ix
=
1
;
ix
<
m_nelx
;
++
ix
)
out
(
ix
-
1
)
=
ix
+
m_nely
*
(
m_nelx
+
1
);
return
out
;
}
// ----------------------- node-numbers along the left edge, without corners -----------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesLeftOpenEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nely
-
1
});
for
(
size_t
iy
=
1
;
iy
<
m_nely
;
++
iy
)
out
(
iy
-
1
)
=
iy
*
(
m_nelx
+
1
);
return
out
;
}
// ---------------------- node-numbers along the right edge, without corners -----------------------
inline
xt
::
xtensor
<
size_t
,
1
>
Regular
::
nodesRightOpenEdge
()
const
{
xt
::
xtensor
<
size_t
,
1
>
out
=
xt
::
empty
<
size_t
>
({
m_nely
-
1
});
for
(
size_t
iy
=
1
;
iy
<
m_nely
;
++
iy
)
out
(
iy
-
1
)
=
iy
*
(
m_nelx
+
1
)
+
m_nelx
;
return
out
;
}
// ----------------------------- node-number of the bottom-left corner -----------------------------
inline
size_t
Regular
::
nodesBottomLeftCorner
()
const
{
return
0
;
}
// ---------------------------- node-number of the bottom-right corner -----------------------------
inline
size_t
Regular
::
nodesBottomRightCorner
()
const
{
return
m_nelx
;
}
// ------------------------------ node-number of the top-left corner -------------------------------
inline
size_t
Regular
::
nodesTopLeftCorner
()
const
{
return
m_nely
*
(
m_nelx
+
1
);
}
// ------------------------------ node-number of the top-right corner ------------------------------
inline
size_t
Regular
::
nodesTopRightCorner
()
const
{
return
m_nely
*
(
m_nelx
+
1
)
+
m_nelx
;
}
// ----------------------------- node-number of the corners (aliases) ------------------------------
inline
size_t
Regular
::
nodesLeftBottomCorner
()
const
{
return
nodesBottomLeftCorner
();
}
inline
size_t
Regular
::
nodesLeftTopCorner
()
const
{
return
nodesTopLeftCorner
();
}
inline
size_t
Regular
::
nodesRightBottomCorner
()
const
{
return
nodesBottomRightCorner
();
}
inline
size_t
Regular
::
nodesRightTopCorner
()
const
{
return
nodesTopRightCorner
();
}
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline
xt
::
xtensor
<
size_t
,
2
>
Regular
::
nodesPeriodic
()
const
{
// edges (without corners)
xt
::
xtensor
<
size_t
,
1
>
bot
=
nodesBottomOpenEdge
();
xt
::
xtensor
<
size_t
,
1
>
top
=
nodesTopOpenEdge
();
xt
::
xtensor
<
size_t
,
1
>
lft
=
nodesLeftOpenEdge
();
xt
::
xtensor
<
size_t
,
1
>
rgt
=
nodesRightOpenEdge
();
// allocate nodal ties
// - number of tying per category
size_t
tedge
=
bot
.
size
()
+
lft
.
size
();
size_t
tnode
=
3
;
// - allocate
xt
::
xtensor
<
size_t
,
2
>
out
=
xt
::
empty
<
size_t
>
({
tedge
+
tnode
,
std
::
size_t
(
2
)});
// counter
size_t
i
=
0
;
// tie all corners to one corner
out
(
i
,
0
)
=
nodesBottomLeftCorner
();
out
(
i
,
1
)
=
nodesBottomRightCorner
();
++
i
;
out
(
i
,
0
)
=
nodesBottomLeftCorner
();
out
(
i
,
1
)
=
nodesTopRightCorner
();
++
i
;
out
(
i
,
0
)
=
nodesBottomLeftCorner
();
out
(
i
,
1
)
=
nodesTopLeftCorner
();
++
i
;
// tie all corresponding edges to each other
for
(
size_t
j
=
0
;
j
<
bot
.
size
()
;
++
j
){
out
(
i
,
0
)
=
bot
(
j
);
out
(
i
,
1
)
=
top
(
j
);
++
i
;
}
for
(
size_t
j
=
0
;
j
<
lft
.
size
()
;
++
j
){
out
(
i
,
0
)
=
lft
(
j
);
out
(
i
,
1
)
=
rgt
(
j
);
++
i
;
}
return
out
;
}
// ------------------------------ node-number that lies in the origin ------------------------------
inline
size_t
Regular
::
nodesOrigin
()
const
{
return
nodesBottomLeftCorner
();
}
// ------------------------- DOF numbers per node (sequentially numbered) --------------------------
inline
xt
::
xtensor
<
size_t
,
2
>
Regular
::
dofs
()
const
{
return
GooseFEM
::
Mesh
::
dofs
(
m_nnode
,
m_ndim
);
}
// ------------------------ DOP-numbers with periodic dependencies removed -------------------------
inline
xt
::
xtensor
<
size_t
,
2
>
Regular
::
dofsPeriodic
()
const
{
// DOF-numbers for each component of each node (sequential)
xt
::
xtensor
<
size_t
,
2
>
out
=
GooseFEM
::
Mesh
::
dofs
(
m_nnode
,
m_ndim
);
// periodic node-pairs
xt
::
xtensor
<
size_t
,
2
>
nodePer
=
nodesPeriodic
();
// eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs
for
(
size_t
i
=
0
;
i
<
nodePer
.
shape
()[
0
]
;
++
i
)
for
(
size_t
j
=
0
;
j
<
m_ndim
;
++
j
)
out
(
nodePer
(
i
,
1
),
j
)
=
out
(
nodePer
(
i
,
0
),
j
);
// renumber "out" to be sequential
return
GooseFEM
::
Mesh
::
renumber
(
out
);
}
// ------------------------------ get the orientation of each element ------------------------------
inline
xt
::
xtensor
<
int
,
1
>
getOrientation
(
const
xt
::
xtensor
<
double
,
2
>
&
coor
,
const
xt
::
xtensor
<
size_t
,
2
>
&
conn
)
{
assert
(
conn
.
shape
()[
1
]
==
3
);
assert
(
coor
.
shape
()[
1
]
==
2
);
double
k
;
size_t
nelem
=
conn
.
shape
()[
0
];
xt
::
xtensor
<
int
,
1
>
out
=
xt
::
empty
<
int
>
({
nelem
});
for
(
size_t
ielem
=
0
;
ielem
<
nelem
;
++
ielem
)
{
auto
v1
=
xt
::
view
(
coor
,
conn
(
ielem
,
0
),
xt
::
all
())
-
xt
::
view
(
coor
,
conn
(
ielem
,
1
),
xt
::
all
());
auto
v2
=
xt
::
view
(
coor
,
conn
(
ielem
,
2
),
xt
::
all
())
-
xt
::
view
(
coor
,
conn
(
ielem
,
1
),
xt
::
all
());
k
=
v1
(
0
)
*
v2
(
1
)
-
v2
(
0
)
*
v1
(
1
);
if
(
k
<
0
)
out
(
ielem
)
=
-
1
;
else
out
(
ielem
)
=
+
1
;
}
return
out
;
}
// ------------------------------ set the orientation of each element ------------------------------
inline
xt
::
xtensor
<
size_t
,
2
>
setOrientation
(
const
xt
::
xtensor
<
double
,
2
>
&
coor
,
const
xt
::
xtensor
<
size_t
,
2
>
&
conn
,
int
orientation
)
{
assert
(
conn
.
shape
()[
1
]
==
3
);
assert
(
coor
.
shape
()[
1
]
==
2
);
assert
(
orientation
==
-
1
||
orientation
==
+
1
);
xt
::
xtensor
<
int
,
1
>
val
=
getOrientation
(
coor
,
conn
);
return
setOrientation
(
coor
,
conn
,
val
,
orientation
);
}
// -------------------- set the orientation of each element to a certain value ---------------------
inline
xt
::
xtensor
<
size_t
,
2
>
setOrientation
(
const
xt
::
xtensor
<
double
,
2
>
&
coor
,
const
xt
::
xtensor
<
size_t
,
2
>
&
conn
,
const
xt
::
xtensor
<
int
,
1
>
&
val
,
int
orientation
)
{
assert
(
conn
.
shape
()[
1
]
==
3
);
assert
(
coor
.
shape
()[
1
]
==
2
);
assert
(
conn
.
shape
()[
0
]
==
val
.
size
()
);
assert
(
orientation
==
-
1
||
orientation
==
+
1
);
// avoid compiler warning
UNUSED
(
coor
);
size_t
nelem
=
conn
.
shape
()[
0
];
xt
::
xtensor
<
size_t
,
2
>
out
=
conn
;
for
(
size_t
ielem
=
0
;
ielem
<
nelem
;
++
ielem
)
if
(
(
orientation
==
-
1
and
val
(
ielem
)
>
0
)
or
(
orientation
==
+
1
and
val
(
ielem
)
<
0
)
)
std
::
swap
(
out
(
ielem
,
2
)
,
out
(
ielem
,
1
)
);
return
out
;
}
// ------------------------------------- re-triangulation API --------------------------------------
inline
xt
::
xtensor
<
size_t
,
2
>
retriangulate
(
const
xt
::
xtensor
<
double
,
2
>
&
coor
,
const
xt
::
xtensor
<
size_t
,
2
>
&
conn
,
int
orientation
)
{
// get the orientation of all elements
xt
::
xtensor
<
int
,
1
>
dir
=
getOrientation
(
coor
,
conn
);
// check the orientation
bool
eq
=
static_cast
<
size_t
>
(
std
::
abs
(
xt
::
sum
(
dir
)[
0
]))
==
conn
.
shape
()[
0
];
// new connectivity
xt
::
xtensor
<
size_t
,
2
>
out
;
// perform re-triangulation
// - use "TriUpdate"
if
(
eq
)
{
Private
::
TriUpdate
obj
(
coor
,
conn
);
obj
.
eval
();
out
=
obj
.
conn
();
}
// - using TriCompute
else
{
throw
std
::
runtime_error
(
"Work-in-progress, has to be re-triangulated using 'TriCompute'"
);
}
return
setOrientation
(
coor
,
out
,
orientation
);
}
// ================================= GooseFEM::Mesh::Tri3::Private =================================
namespace
Private
{
// ------------------------------------------ constructor ------------------------------------------
inline
TriUpdate
::
TriUpdate
(
const
xt
::
xtensor
<
double
,
2
>
&
coor
,
const
xt
::
xtensor
<
size_t
,
2
>
&
conn
)
:
m_conn
(
conn
),
m_coor
(
coor
)
{
assert
(
conn
.
shape
()[
1
]
==
3
);
assert
(
coor
.
shape
()[
1
]
==
2
);
// store shapes
m_nnode
=
coor
.
shape
()[
0
];
m_ndim
=
coor
.
shape
()[
1
];
m_nelem
=
conn
.
shape
()[
0
];
m_nne
=
conn
.
shape
()[
1
];
// set default to out-of-bounds, to make clear that nothing happened yet
m_elem
=
m_nelem
*
xt
::
ones
<
size_t
>
({
2
});
m_node
=
m_nnode
*
xt
::
ones
<
size_t
>
({
4
});
edge
();
}
// -------------------------- compute neighbors per edge of all elements ---------------------------
inline
void
TriUpdate
::
edge
()
{
// signal that nothing has been set
m_edge
=
m_nelem
*
xt
::
ones
<
size_t
>
({
m_nelem
,
m_nne
});
std
::
vector
<
size_t
>
idx
=
{
0
,
1
,
2
};
// lists to convert connectivity -> edges
std
::
vector
<
size_t
>
jdx
=
{
1
,
2
,
0
};
// lists to convert connectivity -> edges
std
::
vector
<
Edge
>
edge
;
edge
.
reserve
(
m_nelem
*
idx
.
size
());
for
(
size_t
e
=
0
;
e
<
m_nelem
;
++
e
)
for
(
size_t
i
=
0
;
i
<
m_nne
;
++
i
)
edge
.
push_back
(
Edge
(
m_conn
(
e
,
idx
[
i
]),
m_conn
(
e
,
jdx
[
i
])
,
e
,
i
,
true
)
);
std
::
sort
(
edge
.
begin
()
,
edge
.
end
()
,
Edge_sort
);
for
(
size_t
i
=
0
;
i
<
edge
.
size
()
-
1
;
++
i
)
{
if
(
edge
[
i
].
n1
==
edge
[
i
+
1
].
n1
and
edge
[
i
].
n2
==
edge
[
i
+
1
].
n2
)
{
m_edge
(
edge
[
i
].
elem
,
edge
[
i
].
edge
)
=
edge
[
i
+
1
].
elem
;
m_edge
(
edge
[
i
+
1
].
elem
,
edge
[
i
+
1
].
edge
)
=
edge
[
i
].
elem
;
}
}
}
// ---------------------------- update edges around renumbered elements ----------------------------
inline
void
TriUpdate
::
chedge
(
size_t
edge
,
size_t
old_elem
,
size_t
new_elem
)
{
size_t
m
;
size_t
neigh
=
m_edge
(
old_elem
,
edge
);
if
(
neigh
>=
m_nelem
)
return
;
for
(
m
=
0
;
m
<
m_nne
;
++
m
)
if
(
m_edge
(
neigh
,
m
)
==
old_elem
)
break
;
m_edge
(
neigh
,
m
)
=
new_elem
;
}
// --------------------------------- re-triangulate the full mesh ----------------------------------
inline
bool
TriUpdate
::
eval
()
{
bool
change
=
false
;
while
(
increment
()
)
{
change
=
true
;
}
return
change
;
}
// ----------------------------------- one re-triangulation step -----------------------------------
inline
bool
TriUpdate
::
increment
()
{
size_t
ielem
,
jelem
,
iedge
,
jedge
;
double
phi1
,
phi2
;
xt
::
xtensor_fixed
<
size_t
,
xt
::
xshape
<
4
>>
c
=
xt
::
empty
<
size_t
>
({
4
});
xt
::
xtensor_fixed
<
size_t
,
xt
::
xshape
<
4
>>
n
=
xt
::
empty
<
size_t
>
({
4
});
// loop over all elements
for
(
ielem
=
0
;
ielem
<
m_nelem
;
++
ielem
)
{
// loop over all edges
for
(
iedge
=
0
;
iedge
<
m_nne
;
++
iedge
)
{
// only proceed if the edge is shared with another element
if
(
m_edge
(
ielem
,
iedge
)
>=
m_nelem
)
continue
;
// read "jelem"
jelem
=
m_edge
(
ielem
,
iedge
);
// find the edge involved for "jelem"
for
(
jedge
=
0
;
jedge
<
m_nne
;
++
jedge
)
if
(
m_edge
(
jelem
,
jedge
)
==
ielem
)
break
;
// convert to four static nodes
// - read first three from "ielem"
if
(
iedge
==
0
)
{
c
(
0
)
=
m_conn
(
ielem
,
2
);
c
(
1
)
=
m_conn
(
ielem
,
0
);
c
(
2
)
=
m_conn
(
ielem
,
1
);
}
else
if
(
iedge
==
1
)
{
c
(
0
)
=
m_conn
(
ielem
,
0
);
c
(
1
)
=
m_conn
(
ielem
,
1
);
c
(
2
)
=
m_conn
(
ielem
,
2
);
}
else
if
(
iedge
==
2
)
{
c
(
0
)
=
m_conn
(
ielem
,
1
);
c
(
1
)
=
m_conn
(
ielem
,
2
);
c
(
2
)
=
m_conn
(
ielem
,
0
);
}
// - read last from "jelem"
if
(
jedge
==
0
)
{
c
(
3
)
=
m_conn
(
jelem
,
2
);
}
else
if
(
jedge
==
1
)
{
c
(
3
)
=
m_conn
(
jelem
,
0
);
}
else
if
(
jedge
==
2
)
{
c
(
3
)
=
m_conn
(
jelem
,
1
);
}
// construct edge vectors
auto
a1
=
xt
::
view
(
m_coor
,
c
(
1
),
xt
::
all
())
-
xt
::
view
(
m_coor
,
c
(
0
),
xt
::
all
());
auto
b1
=
xt
::
view
(
m_coor
,
c
(
2
),
xt
::
all
())
-
xt
::
view
(
m_coor
,
c
(
0
),
xt
::
all
());
auto
a2
=
xt
::
view
(
m_coor
,
c
(
1
),
xt
::
all
())
-
xt
::
view
(
m_coor
,
c
(
3
),
xt
::
all
());
auto
b2
=
xt
::
view
(
m_coor
,
c
(
2
),
xt
::
all
())
-
xt
::
view
(
m_coor
,
c
(
3
),
xt
::
all
());
// compute angles of the relevant corners
phi1
=
std
::
acos
((
a1
(
0
)
*
b1
(
0
)
+
a1
(
1
)
*
b1
(
1
))
/
(
std
::
pow
((
a1
(
0
)
*
a1
(
0
)
+
a1
(
1
)
*
a1
(
1
)),
0.5
)
*
std
::
pow
((
b1
(
0
)
*
b1
(
0
)
+
b1
(
1
)
*
b1
(
1
)),
0.5
)));
phi2
=
std
::
acos
((
a2
(
0
)
*
b2
(
0
)
+
a2
(
1
)
*
b2
(
1
))
/
(
std
::
pow
((
a2
(
0
)
*
a2
(
0
)
+
a2
(
1
)
*
a2
(
1
)),
0.5
)
*
std
::
pow
((
b2
(
0
)
*
b2
(
0
)
+
b2
(
1
)
*
b2
(
1
)),
0.5
)));
// update mesh if needed
if
(
phi1
+
phi2
>
M_PI
)
{
// update connectivity
m_conn
(
ielem
,
0
)
=
c
(
0
);
m_conn
(
ielem
,
1
)
=
c
(
3
);
m_conn
(
ielem
,
2
)
=
c
(
2
);
m_conn
(
jelem
,
0
)
=
c
(
0
);
m_conn
(
jelem
,
1
)
=
c
(
1
);
m_conn
(
jelem
,
2
)
=
c
(
3
);
// change list with neighbors for the elements around (only two neighbors change)
if
(
iedge
==
0
)
{
chedge
(
2
,
ielem
,
jelem
);
}
else
if
(
iedge
==
1
)
{
chedge
(
0
,
ielem
,
jelem
);
}
else
if
(
iedge
==
2
)
{
chedge
(
1
,
ielem
,
jelem
);
}
if
(
jedge
==
0
)
{
chedge
(
2
,
jelem
,
ielem
);
}
else
if
(
jedge
==
1
)
{
chedge
(
0
,
jelem
,
ielem
);
}
else
if
(
jedge
==
2
)
{
chedge
(
1
,
jelem
,
ielem
);
}
// convert to four static nodes
if
(
iedge
==
0
)
{
n
(
0
)
=
m_edge
(
ielem
,
2
);
n
(
3
)
=
m_edge
(
ielem
,
1
);
}
else
if
(
iedge
==
1
)
{
n
(
0
)
=
m_edge
(
ielem
,
0
);
n
(
3
)
=
m_edge
(
ielem
,
2
);
}
else
if
(
iedge
==
2
)
{
n
(
0
)
=
m_edge
(
ielem
,
1
);
n
(
3
)
=
m_edge
(
ielem
,
0
);
}
if
(
jedge
==
0
)
{
n
(
1
)
=
m_edge
(
jelem
,
1
);
n
(
2
)
=
m_edge
(
jelem
,
2
);
}
else
if
(
jedge
==
1
)
{
n
(
1
)
=
m_edge
(
jelem
,
2
);
n
(
2
)
=
m_edge
(
jelem
,
0
);
}
else
if
(
jedge
==
2
)
{
n
(
1
)
=
m_edge
(
jelem
,
0
);
n
(
2
)
=
m_edge
(
jelem
,
1
);
}
// store the neighbors for the changed elements
m_edge
(
ielem
,
0
)
=
jelem
;
m_edge
(
jelem
,
0
)
=
n
(
0
)
;
m_edge
(
ielem
,
1
)
=
n
(
2
)
;
m_edge
(
jelem
,
1
)
=
n
(
1
)
;
m_edge
(
ielem
,
2
)
=
n
(
3
)
;
m_edge
(
jelem
,
2
)
=
ielem
;
// store information for transfer algorithm
m_node
=
c
;
m_elem
(
0
)
=
ielem
;
m_elem
(
1
)
=
jelem
;
return
true
;
}
}
}
return
false
;
}
// ------------------------------------------ constructor ------------------------------------------
inline
Edge
::
Edge
(
size_t
i
,
size_t
j
,
size_t
el
,
size_t
ed
,
bool
sort
)
:
n1
(
i
),
n2
(
j
),
elem
(
el
),
edge
(
ed
)
{
if
(
sort
&&
n1
>
n2
)
std
::
swap
(
n1
,
n2
);
}
// --------------------------------------- compare two edges ---------------------------------------
inline
bool
Edge_cmp
(
Edge
a
,
Edge
b
)
{
if
(
a
.
n1
==
b
.
n1
and
a
.
n2
==
b
.
n2
)
return
true
;
return
false
;
}
// ----------------------- sort edges by comparing the first and second node -----------------------
inline
bool
Edge_sort
(
Edge
a
,
Edge
b
)
{
if
(
a
.
n1
<
b
.
n1
or
a
.
n2
<
b
.
n2
)
return
true
;
return
false
;
}
// -------------------------------------------------------------------------------------------------
}
// namespace Private
// -------------------------------------------------------------------------------------------------
}}}
// namespace ...
// =================================================================================================
#endif
Event Timeline
Log In to Comment