Commit f9b3cffd authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

And moved Reinit library to src and updated makefiles.

parent 1829f405
SHELL = /bin/sh
COMPILE = g++
DEFS = -DPACKAGE=\"AMDiS\" -DVERSION=\"0.1\" -DHAVE_DLFCN_H=1
DEBUG = 0
AMDIS_DIR = ../
INCLUDES = -I./src -I$(AMDIS_DIR)/src -I$(AMDIS_DIR)/compositeFEM/src -I$(AMDIS_DIR)/lib/mtl4
ifeq ($(strip $(DEBUG)), 0)
CPPFLAGS = -O2
else
CPPFLAGS = -g -O0
endif
CPPFLAGS += -ftemplate-depth-40
VPATH = .:$(AMDIS_DIR)/Reinit/src
LIBOFILES = \
BoundaryElementDist.o \
BoundaryElementNormalDist.o \
BoundaryElementLevelSetDist.o \
BoundaryElementTopDist.o \
BoundaryElementEdgeDist.o \
ElementLevelSet.o \
ElementUpdate_2d.o \
ElementUpdate_3d.o \
HL_SignedDist.o \
HL_SignedDistTraverse.o \
VelocityExt.o \
VelocityExtFromVelocityField.o \
NormEps.o
all: libreinit.a
# statische Bibliothek
libreinit.a: $(LIBOFILES)
rm -f lib/libreinit.a
ar cq lib/libreinit.a $(LIBOFILES)
.cc.o: $*.cc
$(COMPILE) $(DEFS) $(INCLUDES) $(CPPFLAGS) -c -o $*.o $<
clean:
-rm -rf *.o
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.
#include "HL_SignedDistLevels.h"
#include "VelocityExtFromVelocityField.h"
void
HL_SignedDistLevels::initializeBoundary()
{
FUNCNAME("HL_SignedDistLevels::initializeBoundary()");
TraverseStack stack;
FixVec<double, VERTEX> distVec(dim, NO_INIT);
int elStatus;
const int *elVertStatusVec;
int NumVertIntPoints;
int ElNum;
clock_t time1;
clock_t time2;
double TIME;
//=== transvering Mesh to SMI and add quantities
time1 = clock();
Mesh_to_SMI_and_quantity();
time2 = clock();
TIME = TIME_USED(time1, time2);
cout <<"Zeit zum Transformieren nach SMI: "<<TIME<<" sec.\n";
// ===== All non-boundary vertices are initialized with "infinity". =====
sD_DOF->set(inftyValue);
// ===== Traverse mesh and initialize boundary elements. =====
const int nBasFcts = feSpace->getBasisFcts()->getNumber();
DegreeOfFreedom *locInd = GET_MEMORY(DegreeOfFreedom, nBasFcts);
ElInfo *elInfo;
if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
elInfo = stack.traverseFirst(feSpace->getMesh(),
-1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_BOUND |
Mesh::FILL_COORDS |
Mesh::FILL_GRD_LAMBDA);
}
else {
elInfo = stack.traverseFirst(feSpace->getMesh(),
-1,
Mesh::CALL_LEAF_EL |
Mesh::FILL_BOUND |
Mesh::FILL_COORDS);
}
while(elInfo) {
// Set elInfo in case velocity extension from velocity field is used.
if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
((VelocityExtFromVelocityField *)velExt)->setElInfo(elInfo);
}
// Get local indices of vertices.
feSpace->getBasisFcts()->getLocalIndices(
const_cast<Element *>(elInfo->getElement()),
const_cast<DOFAdmin *>(feSpace->getAdmin()),
locInd);
// Get element status.
elStatus = elLS->createElementLevelSet(elInfo);
//Get some information for creating the first list
elVertStatusVec = elLS->getElVertStatusVec();
//Beginn creating the first list
NumVertIntPoints = elLS->getNumVertIntPoints();
ElNum = elInfo->getElement()->getIndex();
createFirstList( elStatus, elVertStatusVec, ElNum, NumVertIntPoints);
//end creating the first list
// Is element cut by the interface ?
if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
// Reset element distance vector.
for (int i=0; i<=dim; ++i) {
distVec[i] = inftyValue;
}
// Mark all vertices as boundary vertices.
for (int i=0; i<=dim; ++i) {
(*bound_DOF)[locInd[i]] = 1.0;
}
// Calculate distance for all vertices.
if (bndElDist->calcDistOnBoundaryElement(elInfo, distVec) !=
ElementLevelSet::LEVEL_SET_BOUNDARY) {
ERROR_EXIT("error in distance calculation !\n");
}
else {
// If distance is smaller, correct to new distance.
for (int i=0; i<=dim; ++i) {
if ((*sD_DOF)[locInd[i]] > distVec[i]) {
(*sD_DOF)[locInd[i]] = distVec[i];
//If distance is corrected, calculate new velocity.
if(velExt != NULL)
{
velExt->calcVelocityBoundary(locInd, i);
}
}
}
}
} // end of: elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY
else if (elLS->getNumVertIntPoints() != 0){
// // Interface cuts element only in vertices.
// elVertStatusVec = elLS->getElVertStatusVec();
// for (int i=0; i<=dim; ++i) {
// if (elVertStatusVec[i] == ElementLevelSet::LEVEL_SET_BOUNDARY) {
// (*bound_DOF)[locInd[i]] = 1.0;
// (*sD_DOF)[locInd[i]] = 0.0;
// }
// }
ERROR_EXIT("intersection point is element vertex. should never happen !\n");
}
elInfo = stack.traverseNext(elInfo);
} // end of: mesh traverse
FREE_MEMORY(locInd, DegreeOfFreedom, nBasFcts);
return;
}
void
HL_SignedDistLevels::Mesh_to_SMI_and_quantity()
{
if(!smiAdapter)
{
smiAdapter = NEW SMIAdapter(1, 1,
const_cast<FiniteElemSpace*>(feSpace),
1, -1);
}
// cout << "\n\n\tSMI-Adapter angelegt !\n\n";
//====== transfer Mesh to SMI ======================
smiAdapter->addNeighbourInfo();
// cout << "\n\n\tNachbarschafts-Infos in SMI ergaenzt !\n\n";
smiAdapter-> transferMeshToSMI();
// cout << "\n\n\tGitter nach SMI geschrieben !\n\n";
int smiError = SMI_Begin_write_transaction(1,1);
TEST_EXIT(smiError == SMI_OK)
("SMI_Begin_write_transaction() failed with error %d\n", smiError);
int nul = 0;
int *null = new int [dim+1];
for (int i=0; i<=dim; i++)
{
null[i] = nul;
}
//which pair of element and node is saved in list
smiError =
SMI_Add_quantity(1, 1, 1, SMI_LOC_ELEM, SMI_TYPE_INT, dim+1, null);
TEST_EXIT(smiError == SMI_OK)
("SMI_Add_quantity() failed with error %d\n", smiError);
//which node is a boundary-node
smiError = SMI_Add_quantity(1, 1, 2, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
TEST_EXIT(smiError == SMI_OK)
("SMI_Add_quantity() failed with error %d\n", smiError);
//saves the number of tried updates
smiError = SMI_Add_quantity(1, 1, 3, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
TEST_EXIT(smiError == SMI_OK)
("SMI_Add_quantity() failed with error %d\n", smiError);
//saves the number of updates
smiError = SMI_Add_quantity(1, 1, 4, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
TEST_EXIT(smiError == SMI_OK)
("SMI_Add_quantity() failed with error %d\n", smiError);
//how often a node is saved in the second list
smiError = SMI_Add_quantity(1, 1, 5, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
TEST_EXIT(smiError == SMI_OK)
("SMI_Add_quantity() failed with error %d\n", smiError);
//saves the level of an element
smiError = SMI_Add_quantity(1, 1, 6, SMI_LOC_ELEM, SMI_TYPE_INT, 1,
&inftyValue);
TEST_EXIT(smiError == SMI_OK)
("SMI_Add_quantity() failed with error %d\n", smiError);
// cout << "\n\n\tQuantities in SMI ergaenzt !\n\n";
delete [] null;
}
void
HL_SignedDistLevels::createFirstList(const int elStatus,
const int *elVertStatusVec,
int ElNum,
const int NumVertIntPoints)
{
int *sv = new int[dim+1];
int *nodeIndices;
int nul = 0;
int smiError;
if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY)
{
List_Element.push(ElNum);
smiError =
SMI_Set_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, 1, &ElNum, &nul);
TEST_EXIT(smiError == SMI_OK)
("SMI_Set_quantity_values() failed with error %d\n", smiError);
//saving which nodes are boundary-nodes
for (int i=0; i<=dim; i++)
{
sv[i] = 1;
}
smiError =
SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&ElNum),
NULL, &nodeIndices, NULL, NULL);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_elems() failed with error %d\n", smiError);
smiError =
SMI_Set_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 3, nodeIndices, sv);
TEST_EXIT(smiError == SMI_OK)
("SMI_Set_quantity_values() failed with error %d\n", smiError);
}
//if the elemet isn't a boundary-element, but the interface cuts the FE in two nodes
// else if (NumVertIntPoints == dim)
// {
// List_Element.push( ElNum );
// smiError =
// SMI_Set_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, 1, &ElNum, &nul);
// TEST_EXIT(smiError == SMI_OK)
// ("SMI_Set_quantity_values() failed with error %d\n", smiError);
// //saving which nodes are boundary-nodes
// for (int i=0; i<=dim; i++)
// {
// if(elVertStatusVec[i] == ElementLevelSet::LEVEL_SET_BOUNDARY)
// {
// sv[i] = 1;
// }
// else
// {
// sv[i] = 0;
// }
// smiError =
// SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&ElNum),
// NULL, &nodeIndices, NULL, NULL);
// TEST_EXIT(smiError == SMI_OK)
// ("SMI_Get_elems() failed with error %d\n", smiError);
// smiError =
// SMI_Set_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 3,
// nodeIndices, sv);
// TEST_EXIT(smiError == SMI_OK)
// ("SMI_Set_quantity_values() failed with error %d\n", smiError);
// }
// }
delete [] sv;
}
void
HL_SignedDistLevels::HL_updateIteration()
{
int numNodes;
int *nodeIndices;
int max_q3 = 0;
int max_q4 = 0;
int sum_q3 = 0;
int sum_q4 = 0;
clock_t time1;
clock_t time2;
double TIME;
int control;
int smiError =
SMI_Get_all_nodes(1, 1, const_cast<DegreeOfFreedom*>(&numNodes),
&nodeIndices);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_all_nodes() failed with error %d\n", smiError);
int *values_q3 = new int [numNodes];
int *values_q4 = new int [numNodes];
GET_PARAMETER(0,"SignedDist->count updates", "%d", &count_updates);
time1 = clock();
control = traverseListElement();
time2 = clock();
TIME = TIME_USED(time1, time2);
cout<<"Zeit zum durchlaufen der ersten (Elementen-) Liste: "<<TIME<<" sec.\n\n";
if(control == 1)
{
createLevels ();
}
/////////////////////////////////
//print chosen level in a file
if (print_level == 1)
{
DOFVector<double> *level_DOF = NEW DOFVector<double>(sD_DOF->getFeSpace(),
"level_DOF");
int numElems;
int *elemIndices;
int *NodeIndices;
smiError = SMI_Get_all_elems(1, 1, &numElems, &elemIndices);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_all_elems() failed with error %d\n", smiError);
int *value_q6 = new int [numElems];
smiError =
SMI_Get_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, numElems,
elemIndices, value_q6);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_quantity_values() failed with error %d\n", smiError);
for(int i=0; i<numElems; i++)
{
smiError =
SMI_Get_elems(1, 1, 1, &elemIndices[i], NULL, &NodeIndices,
NULL, NULL);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_elems() failed with error %d\n", smiError);
for(int j=0; j<=dim; j++)
{
(*level_DOF)[NodeIndices[j]] = -1;
}
}
for(int i=0; i<numElems; i++)
{
smiError =
SMI_Get_elems(1, 1, 1, &elemIndices[i], NULL, &NodeIndices,
NULL, NULL);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_elems() failed with error %d\n", smiError);
if (value_q6[i] == chosen_level)
{
for(int j=0; j<dim+1; j++)
{
(*level_DOF)[NodeIndices[j]] = 0;
}
}
}
FileWriter *levelFileWriter =
NEW FileWriter("level->output",
level_DOF->getFeSpace()->getMesh(),
level_DOF);
levelFileWriter->writeFiles (adaptInfo, false);
DELETE levelFileWriter;
}
//end of the part for printing chosen levels
time1 = clock();
traversingListLevel_ ( sD_DOF );
time2 = clock();
TIME = TIME_USED(time1,time2);
cout<<"Zeit zum Durchlaufen der zweiten Liste (Liste von Listen): "<<TIME<<" sec.\n\n";
std::string smiOutFile;
GET_PARAMETER(0, "SignedDist->count updates->output->filename", &smiOutFile);
cout << "count updates Ausgabe-Datei: " << smiOutFile.c_str() << "\n\n";
ofstream out (smiOutFile.c_str());
smiError =
SMI_Get_quantity_values(1, 1, 3, SMI_TYPE_INT, 1, numNodes, nodeIndices,
values_q3);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_quantity_values() failed with error %d\n", smiError);
smiError =
SMI_Get_quantity_values(1, 1, 4, SMI_TYPE_INT, 1, numNodes, nodeIndices,
values_q4);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_quantity_values() failed with error %d\n", smiError);
for (int i=0; i<numNodes; i++)
{
out<<nodeIndices[i]<<" "<<values_q3[i]<<" "<<values_q4[i]<<"\n";
if(max_q3 < values_q3[i])
{
max_q3 = values_q3[i];
}
if(max_q4 < values_q4[i])
{
max_q4 = values_q4[i];
}
sum_q3 = sum_q3 + values_q3[i];
sum_q4 = sum_q4 + values_q4[i];
}
out<<"\n\n maximale Anzahl an versuchten Updates auf einem Knoten: "<<max_q3<<"\n maximale Anzahl an durchgefuehrten Updates auf einem Knoten: "<<max_q4<<"\n";
out<<"\n Summe aller versuchten Updates: "<<sum_q3<<"\n Summe aller durchgefuehrten updates: "<<sum_q4<<"\n";
out<<"Anzahl an Knoten: "<<numNodes<<"\n";
out.close();
smiError = SMI_End_write_transaction(1, 1);
TEST_EXIT(smiError == SMI_OK)
("SMI_End_write_transaction() failed with error %d\n", smiError);
return;
}
int
HL_SignedDistLevels::traverseListElement()
{
int Element;
bool neighboursExist;
bool neighbourInNewListSet = false;
while (List_Element.size() != 0)
{
Element = List_Element.front();
neighboursExist = collectNeighbours_setLevels (Element, 0, &neighbourInNewListSet);
List_Element.pop();
}
//within the function "collectNeighbours_setLevels" the new pairs are saved in the list "level"
//this list will be one entry of the list "Level_"
if (neighbourInNewListSet)
{
Level_.push_back (Level);
}
while (Level.size() != 0)
{
Level.pop();
}
if (!neighbourInNewListSet)
{
return 0;
}
else
{
return 1;
}
}
void
HL_SignedDistLevels::createLevels ()
{
bool neighbourExists = true;
bool stop = false;
Vert_Struct vertStruct;
int i = 0;
bool elementInNextListSet = true;
int value_q2;
int smiError;
while(elementInNextListSet)
{
stop = false;
neighbourExists = false;
elementInNextListSet = false;
while(Level_[i].size() != 0)
{
vertStruct = Level_[i].front();
neighbourExists = collectNeighbours_setLevels (vertStruct.ElNum, i+1, &elementInNextListSet);
smiError =
SMI_Get_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 1,
&vertStruct.VertNum, &value_q2);
TEST_EXIT(smiError == SMI_OK)
("SMI_Get_quantity_values() failed with error %d\n", smiError);
if(value_q2 == 0)
{
helpLevel.push (vertStruct);
}
Level_[i].pop();
}
Level_[i] = helpLevel;
while (helpLevel.size() != 0)
{
helpLevel.pop();
}
if(elementInNextListSet)
{
Level_.push_back (Level);
}
while (Level.size() != 0)
{
Level.pop();
}
i++;
}
}
bool
HL_SignedDistLevels::collectNeighbours_setLevels (const int Element,
const int currentIndex,
bool *elementInNewListSet)
{
Vert_Struct vertStruct;
int *neighbour=new int[dim+1];
int *oppVertices=new int[dim+1];
int *value = new int [dim+1];
int saved_level;
int value_q5;
int *nodeIndices;
int *nodeIndicesOfElem;
int smiError;