Skip to content

Commit b48a83a

Browse files
authored
Merge pull request #36 from dengwirda/dev
0.9.13 updates: support for non-simplex elements, etc
2 parents 55a014c + 3a80cb3 commit b48a83a

File tree

2 files changed

+110
-153
lines changed

2 files changed

+110
-153
lines changed

src/libcpp/iter_mesh/iter_divs_2.inc

Lines changed: 5 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@
3131
*
3232
--------------------------------------------------------
3333
*
34-
* Last updated: 30 April, 2020
34+
* Last updated: 12 Sept., 2020
3535
*
3636
* Copyright 2013-2020
3737
* Darren Engwirda
@@ -79,8 +79,8 @@
7979
iptr_type static constexpr
8080
_last = pred_type::geom_dims+0 ;
8181

82-
iptr_type static constexpr
83-
_DEG_TRIA3 = (iptr_type)+6 ;
82+
// iptr_type static constexpr
83+
// _DEG_TRIA3 = (iptr_type)+6 ;
8484
// iptr_type static constexpr
8585
// _DEG_QUAD4 = (iptr_type)+4 ;
8686

@@ -152,43 +152,8 @@
152152
POINT_tag, _jset) ;
153153

154154
/*--------------------------------- calc. local topo. */
155-
auto _ideg = _iset.count() ;
156-
auto _ierr =
157-
(iptr_type)(_ideg-_DEG_TRIA3) ;
158-
159-
auto _jdeg = _jset.count() ;
160-
auto _jerr =
161-
(iptr_type)(_jdeg-_DEG_TRIA3) ;
162-
163-
// bail-out early if the topological defect would be
164-
// made worse if the zip is done
165-
166-
if (+1 >std::max(_ierr, _jerr))
167-
return ;
168-
169-
if (+0 >std::min(_ierr, _jerr))
170-
return ;
171-
172-
if (+1 > _ierr + _jerr)
173-
return ;
174-
175-
// if more regular topo. is constructed via the edge
176-
// div, make it easier to do so!
177-
178-
real_type _qerr =
179-
(real_type)-1./9. * _ierr +
180-
(real_type)-1./9. * _jerr ;
181-
182-
real_type _lerr =
183-
(real_type)+5./6. * _lmin ;
184-
185-
_qerr /= std::cbrt(_iout) ; // no oscl. wrt. zip
186-
187-
if (+1 < _ierr + _jerr)
188-
{
189-
_qinc = std::min(_qinc, _qerr);
190-
_lmin = std::min(_lmin, _lerr);
191-
}
155+
if (_iset.count() <= +7) return;
156+
if (_jset.count() <= +7) return;
192157

193158
} // if (_lBIG)
194159

src/libcpp/iter_mesh/iter_mesh_2.hpp

Lines changed: 105 additions & 113 deletions
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@
3131
*
3232
--------------------------------------------------------
3333
*
34-
* Last updated: 30 April, 2020
34+
* Last updated: 12 Sept., 2020
3535
*
3636
* Copyright 2013-2020
3737
* Darren Engwirda
@@ -399,7 +399,7 @@
399399
std::min(_0src, *_iter) ;
400400

401401
_msrc += std::pow(
402-
(real_type)1. / *_iter, +5);
402+
(real_type)1. / *_iter, +7);
403403
}
404404
for (auto _iter = _cdst.head(),
405405
_tend = _cdst.tend();
@@ -410,20 +410,20 @@
410410
std::min(_0dst, *_iter) ;
411411

412412
_mdst += std::pow (
413-
(real_type)1. / *_iter, +5);
413+
(real_type)1. / *_iter, +7);
414414
}
415415

416416
_qtol *= std::max(_0src, _zero);
417417

418418
_msrc = std::pow(
419-
_csrc.count() / _msrc, +1./5.) ;
419+
_csrc.count() / _msrc, +1./7.) ;
420420
_mdst = std::pow(
421-
_cdst.count() / _mdst, +1./5.) ;
421+
_cdst.count() / _mdst, +1./7.) ;
422422

423423
_qtol /=
424-
std::pow(_csrc.count(), 1./5.) ;
424+
std::pow(_csrc.count(), 1./7.) ;
425425
_qtol /=
426-
std::pow(_cdst.count(), 1./5.) ;
426+
std::pow(_cdst.count(), 1./7.) ;
427427

428428
/*---------------------------- test move = "okay" */
429429
if (_0dst >= _GOOD)
@@ -663,7 +663,7 @@
663663
{
664664
/*---------------- optimise single node's coordinates */
665665
iptr_type static
666-
constexpr _ITER = (iptr_type) +5 ;
666+
constexpr _ITER = (iptr_type) +4 ;
667667

668668
_move = (iptr_type)-1 ;
669669

@@ -737,7 +737,7 @@
737737
if (_kern == dqdx_kern) // "relax" dQ./dx
738738
{
739739
real_type _BIAS =
740-
(real_type) +2.0 / 3.0 ;
740+
(real_type) +3.0 / 4.0 ;
741741

742742
for (auto _idim =
743743
pred_type::geom_dims; _idim-- != +0; )
@@ -844,7 +844,7 @@
844844
{
845845
/*---------------- optimise single node's coordinates */
846846
iptr_type static
847-
constexpr _ITER = (iptr_type) +5 ;
847+
constexpr _ITER = (iptr_type) +4 ;
848848

849849
__unreferenced(_geom);
850850
__unreferenced(_hfun);
@@ -1545,13 +1545,16 @@
15451545
{
15461546
public :
15471547
/*------------------------ tuple for edge re-ordering */
1548-
iptr_type _edge ;
1548+
iptr_type _inod ;
1549+
iptr_type _jnod ;
15491550
float _cost ;
15501551
public :
15511552
__inline_call sort_pair (
1552-
iptr_type _esrc ,
1553+
iptr_type _isrc ,
1554+
iptr_type _jsrc ,
15531555
float _csrc
1554-
) : _edge(_esrc), _cost(_csrc) {}
1556+
) : _inod(_isrc), _jnod(_jsrc),
1557+
_cost(_csrc) {}
15551558
} ;
15561559
class sort_less
15571560
{
@@ -1598,138 +1601,127 @@
15981601

15991602
_nzip = +0 ; _ndiv = +0 ;
16001603

1601-
for (auto _node =
1602-
_mesh. node().count(); _node-- != 0; )
1603-
{
1604-
/*--------------------- scan nodes and zip//div edges */
1605-
if (_mesh. node(_node).mark () >= 0 &&
1606-
std::abs (
1607-
_mark._node[_node]) > _imrk - 4 )
1604+
// assemble list of edges attached to "recent" nodes
1605+
1606+
for (auto _enum =
1607+
_mesh. edge().count(); _enum-- != 0; )
16081608
{
1609-
_conn.set_count(+0) ;
1610-
_sort.set_count(+0) ;
1611-
/*----------------------- scan edges adj. to node */
1612-
_mesh.connect_1(
1613-
(iptr_type)_node, POINT_tag, _conn) ;
1609+
auto _eptr =&_mesh.edge(_enum) ;
16141610

1615-
for (auto _edge = _conn.head() ;
1616-
_edge != _conn.tend() ;
1617-
++_edge )
1618-
{
1619-
/*----------------------- sort local edges by len */
1620-
auto _enum = _edge-> _cell ;
1621-
auto _eptr =
1622-
_mesh.edge().head() + _enum ;
1611+
auto _inod = _eptr->node(0) ;
1612+
auto _jnod = _eptr->node(1) ;
16231613

1624-
auto _iptr = _mesh.
1625-
node().head()+_eptr->node(0) ;
1626-
auto _jptr = _mesh.
1627-
node().head()+_eptr->node(1) ;
1614+
auto _iptr = _mesh.
1615+
node().head()+_eptr->node(0) ;
1616+
auto _jptr = _mesh.
1617+
node().head()+_eptr->node(1) ;
16281618

1619+
if (_eptr->mark() >= +0 &&
1620+
( std::abs (
1621+
_mark._node[_inod]) > _imrk - 4 ||
1622+
std::abs (
1623+
_mark._node[_jnod]) > _imrk - 4 ))
1624+
{
16291625
float _lsqr =
16301626
(float)pred_type::length_sq (
16311627
& _iptr->pval(0) ,
16321628
& _jptr->pval(0) ) ;
16331629

16341630
_sort.push_tail(
1635-
sort_pair(_enum, _lsqr)) ;
1631+
sort_pair(_inod, _jnod, _lsqr)) ;
16361632
}
1633+
}
16371634

1638-
if (_sort.empty()) continue;
1635+
if (_sort.empty()) return ;
16391636

1640-
/*----------------------- scan local edges by len */
1641-
algorithms::isort(
1642-
_sort.head() ,
1643-
_sort.tend() , sort_less());
1637+
algorithms::qsort( // sort edge list by lsqr
1638+
_sort.head() ,
1639+
_sort.tend() , sort_less());
16441640

1645-
bool_type _move = false ;
1641+
// scan edges longest-to-shortest and try to div any
1642+
// unvisited edges
16461643

1647-
for (auto _iter = _sort.tail();
1648-
_iter != _sort.hend();
1649-
--_iter )
1650-
{
1651-
/*----------------------- try to "div" local edge */
1652-
auto _eadj = _iter->_edge;
1644+
for (auto _iter = _sort.tail();
1645+
_iter != _sort.hend();
1646+
--_iter )
1647+
{
1648+
/*--------------------------- try to "div" local edge */
1649+
iptr_type _eadj, _enod[2] ;
1650+
_enod[0] = _iter->_inod;
1651+
_enod[1] = _iter->_jnod;
16531652

1654-
auto _eptr =
1655-
_mesh. edge().head()+ _eadj;
1653+
if (MARKNODE(_enod[0])>_imrk) continue ;
1654+
if (MARKNODE(_enod[1])>_imrk) continue ;
16561655

1657-
iptr_type _enod[2] ;
1658-
_enod[0] = _eptr->node(0);
1659-
_enod[1] = _eptr->node(1);
1656+
if (MARKNODE(_enod[0]) < +0 ||
1657+
MARKNODE(_enod[1]) < +0 ) continue ;
16601658

1661-
if (MARKNODE(_enod[0])>_imrk) continue ;
1662-
if (MARKNODE(_enod[1])>_imrk) continue ;
1659+
if(!_mesh.find_edge(
1660+
_enod, _eadj) ) continue ;
16631661

1664-
if (MARKNODE(_enod[0]) < +0 ||
1665-
MARKNODE(_enod[1]) < +0 ) continue ;
1662+
if (_opts.div_())
1663+
{
1664+
/*--------------------------- "div" for topo. + score */
1665+
iptr_type _nnew = -1;
1666+
1667+
bool_type _move;
1668+
_div_edge( _geom, _mesh,
1669+
_hfun, _hval, _opts,
1670+
_imrk, _eadj,
1671+
_kern, _move, _nnew,
1672+
_iset, _jset,
1673+
_aset, _bset,
1674+
_qold, _qnew,
1675+
_qtmp, _QLIM) ;
16661676

1667-
if (_opts.div_())
1677+
if (_move)
16681678
{
1669-
/*----------------------- "div" for topo. + score */
1670-
iptr_type _nnew = -1;
1671-
1672-
_div_edge( _geom, _mesh,
1673-
_hfun, _hval, _opts,
1674-
_imrk, _eadj,
1675-
_kern, _move, _nnew,
1676-
_iset, _jset,
1677-
_aset, _bset,
1678-
_qold, _qnew,
1679-
_qtmp, _QLIM) ;
1680-
1681-
if (_move)
1682-
{
1683-
PUSHMARK; _ndiv += +1; break ;
1684-
}
1679+
PUSHMARK; _ndiv += +1;
16851680
}
16861681
}
1682+
}
16871683

1688-
if (_move) continue ;
1684+
// scan edges shortest-to-longest and try to zip any
1685+
// unvisited edges
16891686

1690-
for (auto _iter = _sort.head();
1691-
_iter != _sort.tend();
1692-
++_iter )
1693-
{
1694-
/*----------------------- try to "zip" local edge */
1695-
auto _eadj = _iter->_edge;
1687+
for (auto _iter = _sort.head();
1688+
_iter != _sort.tend();
1689+
++_iter )
1690+
{
1691+
/*--------------------------- try to "zip" local edge */
1692+
iptr_type _eadj, _enod[2] ;
1693+
_enod[0] = _iter->_inod;
1694+
_enod[1] = _iter->_jnod;
16961695

1697-
auto _eptr =
1698-
_mesh. edge().head()+ _eadj;
1696+
if (MARKNODE(_enod[0])>_imrk) continue ;
1697+
if (MARKNODE(_enod[1])>_imrk) continue ;
16991698

1700-
iptr_type _enod[2] ;
1701-
_enod[0] = _eptr->node(0);
1702-
_enod[1] = _eptr->node(1);
1699+
if (MARKNODE(_enod[0]) < +0 ||
1700+
MARKNODE(_enod[1]) < +0 ) continue ;
17031701

1704-
if (MARKNODE(_enod[0])>_imrk) continue ;
1705-
if (MARKNODE(_enod[1])>_imrk) continue ;
1702+
if(!_mesh.find_edge(
1703+
_enod, _eadj) ) continue ;
17061704

1707-
if (MARKNODE(_enod[0]) < +0 ||
1708-
MARKNODE(_enod[1]) < +0 ) continue ;
1705+
if (_opts.zip_())
1706+
{
1707+
/*--------------------------- "zip" for topo. + score */
1708+
iptr_type _nnew = -1;
1709+
1710+
bool_type _move;
1711+
_zip_edge( _geom, _mesh,
1712+
_hfun, _hval, _opts,
1713+
_eadj,
1714+
_kern, _move, _nnew,
1715+
_iset, _jset,
1716+
_aset, _bset, _cset,
1717+
_qold, _qnew,
1718+
_qtmp, _QLIM) ;
17091719

1710-
if (_opts.zip_())
1720+
if (_move)
17111721
{
1712-
/*----------------------- "zip" for topo. + score */
1713-
iptr_type _nnew = -1;
1714-
1715-
_zip_edge( _geom, _mesh,
1716-
_hfun, _hval, _opts,
1717-
_eadj,
1718-
_kern, _move, _nnew,
1719-
_iset, _jset,
1720-
_aset, _bset, _cset,
1721-
_qold, _qnew,
1722-
_qtmp, _QLIM) ;
1723-
1724-
if (_move)
1725-
{
1726-
PUSHMARK; _nzip += +1; break ;
1727-
}
1722+
PUSHMARK; _nzip += +1;
17281723
}
17291724
}
1730-
1731-
if (_move) continue ;
1732-
}
17331725
}
17341726

17351727
for (auto _iter = _mark._node.head() ;
@@ -1990,7 +1982,7 @@
19901982

19911983
real_type _DLIM =
19921984
(real_type)(1. -
1993-
1. * std::pow(1.0-_QLIM, +2)) ;
1985+
.5 * std::pow(1.0-_QLIM, +2)) ;
19941986

19951987
/*------------------------------ 1. CELL GEOM. PASSES */
19961988

0 commit comments

Comments
 (0)