Commit 02fe7a39 authored by Ansgar Burchardt's avatar Ansgar Burchardt

[!134] build tetrahedron refinement rules into library

Merge branch 'build-tetrahedron-refrules-into-library' into 'master'

ref:staging/dune-uggrid The TET_RULESET compile-time option could be used to
load tetrahedron refinement rules from RefRules.data, however this requires to
locate the file at run-time. This was too inconvenient.

Instead generate RefRules.cc from RefRules.data including the refinement rules
and use that.

See merge request [staging/dune-uggrid!134]

  [staging/dune-uggrid!134]: gitlab.dune-project.org/staging/dune-uggrid/merge_requests/134
parents 8352275f 72b22665
......@@ -4,6 +4,11 @@
[!96]: https://gitlab.dune-project.org/staging/dune-uggrid/merge_requests/96
* The `TET_RULESET` compile-time option has been renamed to
`DUNE_UGGRID_TET_RULESET` and is enabled by default.
It is no longer necessary to provide the `RefRules.data` file.
[!134](https://gitlab.dune-project.org/staging/dune-uggrid/merge_requests/134)
# dune-uggrid 2.6.0 (2018-04-03)
* Transfer of element data during load balancing is supported.
......
......@@ -41,6 +41,13 @@ set(UG_ENABLE_SYSTEM_HEAP ON CACHE BOOL "If ON/True then we are using the operat
set(UG_DDD_MAX_MACROBITS "24" CACHE STRING
"Set number of bits of an unsigned int used to store the process number,
the remaining bits are used to store the local entity id")
set(DUNE_UGGRID_TET_RULESET True CACHE BOOL "Use complete rule set for refinement of tetrahedral elements")
if(TET_RULESET)
set(DUNE_UGGRID_TET_RULESET True)
message(DEPRECATION "The TET_RULESET option has been renamed to DUNE_UGGRID_TET_RULESET")
endif()
if(UG_ENABLE_DEBUGGING)
add_definitions("-DDebug")
set(UG_EXTRAFLAGS "${UG_EXTRAFLAGS} -DDebug")
......
......@@ -20,7 +20,7 @@
#cmakedefine TIME_WITH_SYS_TIME 1
/* Define to 1 if UGGrid should use the complete set of green refinement rules for tetrahedra */
#cmakedefine TET_RULESET 1
#cmakedefine DUNE_UGGRID_TET_RULESET 1
/* Define to 1 if rpc/rpc.h is found (needed for xdr). */
#ifndef HAVE_RPC_RPC_H
......
set(SOURCES algebra.cc enrol.cc evm.cc mgio.cc
ugio.cc ugm.cc cw.cc initgm.cc elements.cc
shapes.cc evalproc.cc rm.cc refine.cc
dlmgr.cc gmcheck.cc er.cc mgheapmgr.cc)
dlmgr.cc gmcheck.cc er.cc mgheapmgr.cc rm-write2file.cc )
# put this "template" into distribution-tarball as well
set(EXTRA_DIST dlmgr.t)
set(gminclude_HEADERS cw.h elements.h gm.h pargm.h evm.h
shapes.h ugm.h dlmgr.h algebra.h rm.h refine.h)
shapes.h ugm.h dlmgr.h algebra.h rm.h refine.h rm-write2file.h)
ug_add_dim_libs(ug_gm OBJECT SOURCES ${SOURCES})# OBJECT_DIM_LIBS gg)
dune_add_test(
NAME rm3-tetrahedron-rules-test
SOURCES rm-tetrahedron-rules-test.cc
COMPILE_DEFINITIONS -DUG_DIM_3
LINK_LIBRARIES ugL ugS3 ${DUNE_LIBS}
CMAKE_GUARD DUNE_UGGRID_TET_RULESET
)
# rm3-show
add_executable(rm3-show rm-show.cc)
target_compile_definitions(rm3-show PRIVATE -DUG_DIM_3)
target_link_libraries(rm3-show PRIVATE ugL ugS3 ${DUNE_LIBS})
# rm3-writeRefRules2file
add_executable(rm3-writeRefRules2file rm-writeRefRules2file.cc rm-write2file.h rm-write2file.cc)
target_compile_definitions(rm3-writeRefRules2file PRIVATE -DUG_DIM_3)
target_link_libraries(rm3-writeRefRules2file PRIVATE ugL ugS3 ${DUNE_LIBS})
install(FILES ${gminclude_HEADERS} DESTINATION ${CMAKE_INSTALL_PKGINCLUDEDIR})
This diff is collapsed.
......@@ -627,7 +627,7 @@ static INT CheckEdge (ELEMENT *theElement, EDGE* theEdge, INT i)
theNode = MIDNODE(theEdge);
if (theNode == NULL)
{
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
if (((REFINE(theElement) == RED) && (TAG(theElement) != TETRAHEDRON))
|| ((TAG(theElement) == TETRAHEDRON) && (NSONS(theElement) == 8)))
#else
......
......@@ -164,7 +164,7 @@ using namespace PPIF;
#define MARKED(e) (MARK(e)!=NO_REFINEMENT)
/* green marked elements were NEWGREEN is true are refined without rule */
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
#define NEWGREEN(e) (TAG(e)==HEXAHEDRON || TAG(e)== PRISM || \
TAG(e)==PYRAMID)
#else
......@@ -374,7 +374,7 @@ static int gridadaptl_timer,ident_timer,overlap_timer,gridcons_timer;
static int algebra_timer;
#endif
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
/* determine number of edge from reduced (i.e. restricted to one side) edgepattern */
/* if there are two edges marked for bisection, if not deliver -1. If the edge- */
/* is not reduced (i.e. marked edges lying on more than one side) deliver -2 */
......@@ -1152,7 +1152,7 @@ static INT ComputePatterns (GRID *theGrid)
SETSIDEPATTERN(theElement,0);
for (i=0; i<SIDES_OF_ELEM(theElement); i++)
{
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
if (CORNERS_OF_SIDE(theElement,i)==4)
{
#endif
......@@ -1160,7 +1160,7 @@ static INT ComputePatterns (GRID *theGrid)
if(SIDE_IN_PATTERN(theElement,thePattern,i))
SETSIDEPATTERN(theElement,
SIDEPATTERN(theElement) | 1<<i);
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
}
#endif
}
......@@ -1181,7 +1181,7 @@ static INT ComputePatterns (GRID *theGrid)
}
#ifdef __THREEDIM__
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
/****************************************************************************/
/*
......@@ -1380,7 +1380,7 @@ static INT CorrectElementSidePattern (ELEMENT *theElement, ELEMENT *theNeighbor,
switch (CORNERS_OF_SIDE(theElement,i))
{
case 3 :
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
/* handle case with 2 edges of the side refined */
if (CorrectTetrahedronSidePattern(theElement,i,theNeighbor,j) != GM_OK)
RETURN(GM_ERROR);
......@@ -1582,7 +1582,7 @@ static INT SetElementRules (GRID *theGrid, ELEMENT *firstElement, INT *cnt)
/* choose best tet_red rule according to (*theFullRefRule)() */
if (TAG(theElement)==TETRAHEDRON && MARKCLASS(theElement)==RED_CLASS)
{
#ifndef TET_RULESET
#ifndef DUNE_UGGRID_TET_RULESET
if ((Mark==TET_RED || Mark==TET_RED_0_5 ||
Mark==TET_RED_1_3))
#endif
......@@ -1868,7 +1868,7 @@ static INT BuildGreenClosure (GRID *theGrid)
/* for pyramids, prisms and hexhedra Patterns2Rules returns 0 */
/* for non red elements, because there is no complete rule set */
/* switch to mark COPY, because COPY rule refines no edges */
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
if (DIM==3 && TAG(theElement)!=TETRAHEDRON)
#else
if (DIM==3)
......@@ -1940,7 +1940,7 @@ static INT BuildGreenClosure (GRID *theGrid)
if (NODE_OF_RULE(theNeighbor,MARK(theNeighbor),
EDGES_OF_ELEM(theNeighbor)+j))
{
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
if (TAG(theNeighbor)==TETRAHEDRON)
printf("ERROR: no side nodes for tetrahedra! side=%d\n",j);
#endif
......@@ -2046,7 +2046,7 @@ static int Gather_ElementInfo (DDD::DDDContext&, DDD_OBJ obj, void *Data)
assert(0); \
}
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
#define COMPARE_MACROX(elem0,elem1,macro,print) \
{ \
INT _mark0,_mark1,_pat0,_pat1; \
......@@ -2111,7 +2111,7 @@ static int Scatter_ElementInfo (DDD::DDDContext&, DDD_OBJ obj, void *Data)
COMPARE_MACROX(theElement,theMaster,MARK,PrintDebug)
COMPARE_MACRO(theElement,theMaster,COARSEN,PrintDebug)
COMPARE_MACRO(theElement,theMaster,USED,PrintDebug)
/* #ifndef TET_RULESET */
/* #ifndef DUNE_UGGRID_TET_RULESET */
#ifdef __THREEDIM__
COMPARE_MACRO(theElement,theMaster,SIDEPATTERN,PrintDebug)
#endif
......@@ -2202,7 +2202,7 @@ static int GridClosure (GRID *theGrid)
do
{
#ifdef __THREEDIM__
#if defined(ModelP) && defined(TET_RULESET)
#if defined(ModelP) && defined(DUNE_UGGRID_TET_RULESET)
/* edge pattern is needed consistently in CorrectTetrahedronSidePattern() */
if (!refine_seq)
{
......@@ -2463,7 +2463,7 @@ INT NS_DIM_PREFIX GetSons (const ELEMENT *theElement, ELEMENT *SonList[MAX_SONS]
static INT RestrictElementMark(ELEMENT *theElement)
{
#ifdef __THREEDIM__
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
EDGE *theEdge;
int j,Rule,Pattern;
#endif
......@@ -2487,7 +2487,7 @@ static INT RestrictElementMark(ELEMENT *theElement)
#endif
#ifdef __THREEDIM__
case TETRAHEDRON :
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
if (MARK(theElement)!=RED)
/** \todo Is REFINE always as red rule available? */
SETMARK(theElement,REFINE(theElement));
......@@ -2526,7 +2526,7 @@ static INT RestrictElementMark(ELEMENT *theElement)
#endif
#ifdef __THREEDIM__
case TETRAHEDRON :
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
/* theElement is not marked from outside, */
/* so find a reg. rule being consistent */
/* with those neighbors of all sons of */
......@@ -3185,7 +3185,7 @@ static int UpdateContext (GRID *theGrid, ELEMENT *theElement, NODE **theElementC
for (INT i=0; i<SIDES_OF_ELEM(theElement); i++)
{
/* no side nodes for triangular sides yet */
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
if (CORNERS_OF_SIDE(theElement,i) == 3) continue;
#endif
......
#include "config.h"
#include <cerrno>
#include <cstring>
#include <dune/common/exceptions.hh>
#include "low/namespace.h"
#include "low/initlow.h"
#include "dev/ugdevices.h"
#include "gm/initgm.h"
#include "gm/rm.h"
#include "gm/rm-write2file.h"
int main(int argc, char** argv)
{
UG::InitLow();
UG::InitDevices();
NS_DIM_PREFIX InitGm();
std::vector<NS_DIM_PREFIX REFRULE> rules;
std::vector<NS_PREFIX SHORT> patterns;
// use an arbitrary additional argument to use the rules from "lib/ugdata/RefRules.data"
if( argc > 1 )
{
// read the file "lib/ugdata/RefRules.data" and use those rules
std::FILE* stream = std::fopen(argv[1], "r");
if (!stream)
DUNE_THROW(Dune::Exception, "Could not open file " << argv[1] << ": " << std::strerror(errno));
NS_DIM_PREFIX readTetrahedronRules(stream, rules, patterns);
std::fclose(stream);
NS_DIM_PREFIX RefRules[NS_DIM_PREFIX TETRAHEDRON] = rules.data();
NS_DIM_PREFIX MaxRules[NS_DIM_PREFIX TETRAHEDRON] = rules.size();
}
// write rules
for(int i=0; i<NS_DIM_PREFIX MaxRules[NS_DIM_PREFIX TETRAHEDRON]; i++)
{
NS_DIM_PREFIX ShowRefRuleX(NS_DIM_PREFIX TETRAHEDRON, i, std::printf);
}
return 0;
}
#include "config.h"
#include <cstdio>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#include "evm.h"
#include "gm.h"
#include "initug.h"
#include "rm.h"
USING_UGDIM_NAMESPACE
USING_UG_NAMESPACE
static bool CheckVolumes(REFRULE const *Rule)
{
bool pass = true;
DOUBLE_VECTOR coords [MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM];
DOUBLE_VECTOR coord;
DOUBLE_VECTOR p[MAX_CORNERS_OF_ELEM];
/* this function work only for TETRAHEDRA!!! */
int tag = TETRAHEDRON;
/* reset coords */
for (int i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
V3_CLEAR(coords[i])
/* compute coordinates of corners */
int j = 0;
for (int i=0; i<CORNERS_OF_REF(tag); i++,j++)
V3_COPY(LOCAL_COORD_OF_REF(tag,i),coords[j])
/* compute coordinates of midnodes */
for (int i=0; i<EDGES_OF_REF(tag); i++,j++)
{
V3_CLEAR(coord)
for (int k=0; k<CORNERS_OF_EDGE; k++)
{
V3_ADD(LOCAL_COORD_OF_REF(tag,CORNER_OF_EDGE_REF(tag,i,k)),coord,coord)
}
V3_SCALE(1.0/CORNERS_OF_EDGE,coord)
V3_COPY(coord,coords[j]);
}
/* compute coordinates of sidenodes */
for (int i=0; i<SIDES_OF_REF(tag); i++,j++)
{
V3_CLEAR(coord)
for (int k=0; k<CORNERS_OF_SIDE_REF(tag,i); k++)
V3_ADD(LOCAL_COORD_OF_REF(tag,CORNER_OF_SIDE_REF(tag,i,k)),coord,coord)
V3_SCALE(1.0/CORNERS_OF_SIDE_REF(tag,i),coord)
V3_COPY(coord,coords[j]);
}
/* compute coordinates of center node */
V3_CLEAR(coord)
for (int i=0; i<CORNERS_OF_REF(tag); i++)
{
V3_ADD(LOCAL_COORD_OF_REF(tag,i),coord,coord)
}
V3_SCALE(1.0/CORNERS_OF_REF(tag),coord)
V3_COPY(coord,coords[j]);
for (int i=0; i<MAX_CORNERS_OF_ELEM+MAX_NEW_CORNERS_DIM; i++)
{
printf("CheckVolumes(): i=%d x=%.8f y=%.8f z=%.8f\n",
i,coords[i][0],coords[i][1],coords[i][2]);
}
/* check the son volumes being really greater than zero */
double sum=0.0;
for (int i=0; i<NSONS_OF_RULE(Rule); i++)
{
for (int k=1; k<4; k++)
{
V3_SUBTRACT(coords[SON_CORNER_OF_RULE(Rule,i,k)],coords[SON_CORNER_OF_RULE(Rule,i,0)],p[k])
}
double sp;
// triple product of son
V3_VECTOR_PRODUCT(p[1],p[2],p[0])
V3_SCALAR_PRODUCT(p[0],p[3],sp)
if (sp<=0.0)
{
printf("negative volume=%f for son=%d rule=%d\n",sp,i,MARK_OF_RULE(Rule));
pass = false;
}
sum += sp;
}
// check if volume of sons add up to volume of original tetrahedra
// 6*Volume_original == sum over 6*volume_sons
// exception for first rule
if( sum != 1.0 && Rule->mark != 0)
{
printf("Rule %d :sum over sons = %f != 1\n",Rule->mark,sum);
pass = false;
}
return pass;
}
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
InitUg(&argc, &argv);
bool pass = true;
std::printf("Testing %d refinement rules for the tetrahedron...\n", MaxRules[TETRAHEDRON]);
REFRULE* rules = RefRules[TETRAHEDRON];
for (int i = 0; i < MaxRules[TETRAHEDRON]; ++i)
pass = pass && CheckVolumes(rules + i);
ExitUg();
return pass ? 0 : 1;
}
This diff is collapsed.
#ifndef RM_WRITE2FILE_H
#define RM_WRITE2FILE_H
// standard C library
#include <cassert>
#include <cstdio>
#include <cmath>
#include <cstdlib>
// std library
#include <stdio.h>
#include <string>
#include <memory>
#include <utility>
#include <iostream>
#include <algorithm>
#include <iterator>
#include <dev/ugdevices.h>
#include "gm/gm.h"
#include "gm/rm.h"
#include "gm/evm.h"
#define TET_RULE_FATHER_SIDE_OFFSET 20
int WriteSonData(std::FILE* stream, NS_DIM_PREFIX sondata const& son);
void WriteRule2File(std::FILE* stream, NS_DIM_PREFIX REFRULE const& rules);
void Write2File(std::FILE* stream, std::vector<NS_DIM_PREFIX REFRULE> const& rules, std::vector<NS_PREFIX SHORT> const& patterns);
START_UGDIM_NAMESPACE
void readTetrahedronRules(FILE* stream, std::vector<NS_DIM_PREFIX REFRULE>& rules, std::vector<NS_PREFIX SHORT>& patterns);
END_UGDIM_NAMESPACE
#endif
#include "config.h"
#include <cerrno>
#include <cstring>
#include <vector>
#include <dune/common/exceptions.hh>
#include "low/namespace.h"
#include "gm/rm-write2file.h"
#include <iostream>
#include <memory>
#include <utility>
int main(int argc, char** argv)
{
if (argc != 3) {
std::cerr << "usage: " << argv[0] << " <RefRules.data> <RefRules.cc>\n";
return 1;
}
std::vector<NS_DIM_PREFIX REFRULE> rules;
std::vector<NS_PREFIX SHORT> patterns;
// Read rules
{
const char* filename = argv[1];
std::FILE* stream = std::fopen(filename, "r");
if (!stream)
DUNE_THROW(Dune::Exception, "Could not open " << filename << " for reading: " << std::strerror(errno));
readTetrahedronRules(stream, rules, patterns);
std::fclose(stream);
}
// Write rules
{
const char* filename = argv[2];
std::FILE* stream = std::fopen(filename, "w");
if (!stream)
DUNE_THROW(Dune::Exception, "Could not open " << filename << " for writing: " << std::strerror(errno));
Write2File(stream, rules, patterns);
if (std::fclose(stream))
DUNE_THROW(Dune::Exception, "E: Closing " << filename << " failed: " << std::strerror(errno));
}
return 0;
}
This diff is collapsed.
......@@ -52,7 +52,7 @@ START_UGDIM_NAMESPACE
/* */
/****************************************************************************/
/* Declaration of TET_RULESET has been moved to config.h */
/* Declaration of DUNE_UGGRID_TET_RULESET has been moved to config.h */
/* uncomment this if you want to use the full rule set for tetrahedra */
/** \brief Defines for edge types */
......@@ -172,7 +172,7 @@ enum {Q_NOREF,
Q_CLOSE_3_3};
#define TET_COPY 1
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_TET_RULESET
#define FULL_REFRULE Pattern2Rule[TETRAHEDRON][0x3F]
#define FULL_REFRULE_0_5 (Pattern2Rule[TETRAHEDRON][0x3F]+1)
#define FULL_REFRULE_1_3 (Pattern2Rule[TETRAHEDRON][0x3F]+2)
......@@ -324,7 +324,7 @@ extern INT MaxNewCorners[TAGS];
extern INT MaxNewEdges[TAGS];
extern INT CenterNodeIndex[TAGS];
extern REFRULE *RefRules[TAGS];
extern std::unique_ptr<SHORT[]> Pattern2Rule[TAGS];
extern const SHORT* Pattern2Rule[TAGS];
extern FULLREFRULEPTR theFullRefRule;
......
......@@ -1304,7 +1304,7 @@ INT NS_DIM_PREFIX GetSideIDFromScratch (ELEMENT *theElement, NODE *theNode)
for (i=0; i<SIDES_OF_ELEM(theFather); i++)
{
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_DUNE_UGGRID_TET_RULESET
if (3 == CORNERS_OF_SIDE(theFather,i)) continue;
#endif
......@@ -1734,7 +1734,7 @@ INT NS_DIM_PREFIX GetNodeContext (const ELEMENT *theElement, NODE **theElementCo
EDGES_OF_ELEM(theElement);
for (i=0; i<SIDES_OF_ELEM(theElement); i++)
{
#ifdef TET_RULESET
#ifdef DUNE_UGGRID_DUNE_UGGRID_TET_RULESET
/* no side nodes for triangular sides yet */
if (CORNERS_OF_SIDE(theElement,i) == 3) continue;
#endif
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment