CompositeFEMMethods.cc 3.89 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// 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.


13
14
15
16
17
18
19
#include "BasisFunction.h"
#include "DOFAdmin.h"
#include "Element.h"
#include "ElInfo.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
20
#include "Initfile.h"
21
22
23
24
25
#include "Traverse.h"

#include "CompositeFEMMethods.h"


26
27
28
void CompositeFEMMethods::setPosLsToVal(DOFVector<double> *dof,
					const double &val,
					const DOFVector<double> *lsFct_dof)
29
30
31
32
33
34
35
36
37
38
{
  DOFVector<double>::Iterator it_dof(dof, USED_DOFS);
  DOFVector<double>::Iterator it_lsFct_dof(
                     const_cast<DOFVector<double> *>(lsFct_dof), 
		     USED_DOFS);
  for (it_dof.reset(), it_lsFct_dof.reset();
       !it_dof.end();
       ++it_dof, ++it_lsFct_dof) {

    // Is vertex in domain with positive level set function values ?
39
    if (*it_lsFct_dof > 0)
40
41
42
43
      *it_dof = val;
  }
}

44
45
46
void CompositeFEMMethods::setPosLsToFct(DOFVector<double> *dof,
					const AbstractFunction<double, WorldVector<double> > *fct,
					const DOFVector<double> *lsFct_dof)
47
{
48
49
50
  const BasisFunction *basisFcts = dof->getFeSpace()->getBasisFcts();
  const DOFAdmin *admin = dof->getFeSpace()->getAdmin();
  const int dim = dof->getFeSpace()->getMesh()->getDim();
51
52
  const DegreeOfFreedom *locInd;

53
54
55
  TraverseStack stack;
  ElInfo *elInfo = stack.traverseFirst(dof->getFeSpace()->getMesh(), -1, 
				       Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
56
57
  ElementVector locVec(basisFcts->getNumber());

58
  while (elInfo) {
59
    const Element *el = elInfo->getElement();
60
61
    
    // Get level set function values for all vertices of element.
62
    lsFct_dof->getLocalVector(el, locVec);
63
64
65
66

    // Get dof indices of vertices.
    locInd = basisFcts->getLocalIndices(el, admin, NULL);

67
    for (int i = 0; i <= dim; i++) {
68
      // Is vertex in domain with positive level set function values ?
69
      if (locVec[i] > 0)
70
71
72
73
74
75
76
	(*dof)[locInd[i]] = (*fct)(elInfo->getCoord(i));
    }

    elInfo = stack.traverseNext(elInfo);
  }
}

77
78
79
void CompositeFEMMethods::printBoundaryElements(const std::string fn_str,
						ElementLevelSet *elLS,
						FiniteElemSpace *feSpace)
80
81
82
83
84
85
{
  FUNCNAME("CompositeFEMMethods::printBoundaryElements()");

  int dim = feSpace->getMesh()->getDim();
  std::string fn_main;
  std::string fn;
86
  Parameters::get(fn_str + "->filename", fn_main);
87
88
89
90
91
92
93
94
95
96
97
  fn = fn_main + "boundary_elements";
  int elStatus;

  ofstream boundaryOut(fn.c_str());

  // ===== Traverse mesh and print boundary elements. =====
  TraverseStack stack;
  int boundEl_cntr = 0;
  WorldVector<double> coord;

  const int nBasFcts = feSpace->getBasisFcts()->getNumber();
98
  DegreeOfFreedom *locInd = new DegreeOfFreedom[nBasFcts];
99
100
101
102
103
104

  ElInfo *loc_elInfo = stack.traverseFirst(feSpace->getMesh(),
					   -1, 
					   Mesh::CALL_LEAF_EL | 
					   Mesh::FILL_BOUND |
					   Mesh::FILL_COORDS);
105
  while (loc_elInfo) {
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

    // Get local indices of vertices.
    feSpace->getBasisFcts()->getLocalIndices(
	  const_cast<Element *>(loc_elInfo->getElement()),
	  const_cast<DOFAdmin *>(feSpace->getAdmin()),
	  locInd); 

    // Get element status.
    elStatus = elLS->createElementLevelSet(loc_elInfo);

    // Is element cut by the interface ?
    if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
      ++boundEl_cntr;

      boundaryOut << loc_elInfo->getElement()->getIndex() << ": \t";

      for (int i=0; i<=dim; ++i) {

	coord = loc_elInfo->getCoord(i);

	if (i > 0) {
	  boundaryOut << "\n\t";
	}

	boundaryOut << "\t" << coord[0];
	for (int j=1; j<dim; ++j) {
	  boundaryOut << "; " << coord[j] ;
	}
	boundaryOut << " (" << locInd[i] << ")";
      }
      boundaryOut << "\n";
    }

    loc_elInfo = stack.traverseNext(loc_elInfo);

  }  // end of: mesh traverse

143
  delete [] locInd;
144
145
146
147
148

  boundaryOut << "\nNumber of boundary elements: \t" << boundEl_cntr << "\n";

  boundaryOut.close();
}