Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

Helpers.cc 14.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


19 20 21 22
#include "Helpers.h"

using namespace AMDiS;

Praetorius, Simon's avatar
Helpers  
Praetorius, Simon committed
23
namespace Helpers {
24

25
  void calcMaxOnXAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima)
Praetorius, Simon's avatar
Praetorius, Simon committed
26
  {
27 28
    std::vector<WorldVector<double> > x;
    std::vector<double> y;
29

30 31 32
    double dimX=1.0,dimY=1.0;
    Parameters::get("mesh->dimension[x]",dimX);
    Parameters::get("mesh->dimension[y]",dimY);
33

34 35 36
    interpolOverLine(*rho,-dimX,0.0,dimX,0.0,1000,x,y);
    if (x.size() != y.size())
      throw(std::runtime_error("x.size /= y.size! Sollte nicht passieren!"));
37

38 39 40 41 42
    if(x.size()>0) {
      for(unsigned i=1; i<x.size()-1; ++i) {
	if(y[i-1]<y[i] && y[i+1]<y[i]) {
	  maxima.push_back(std::make_pair(x[i],y[i]));
	}
43 44 45 46
      }
    }
  }

47 48
  
  void calcMaxOnYAxis(DOFVector<double> *rho, std::vector<std::pair<WorldVector<double>, double> > &maxima)
Praetorius, Simon's avatar
Praetorius, Simon committed
49
  {
50

51 52
    std::vector<WorldVector<double> > x;
    std::vector<double> y;
53

54 55 56
    double dimX,dimY;
    Parameters::get("mesh->dimension[x]",dimX);
    Parameters::get("mesh->dimension[y]",dimY);
57

58
    interpolOverLine(*rho,0.0,-dimY,0.0,dimY,1000,x,y);
59

60 61 62 63 64 65 66 67 68
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (MPI::COMM_WORLD.Get_rank() == 0)
    #endif
    {
    for(unsigned i=1; i<x.size()-1; ++i) {
      if(y[i-1]<y[i] && y[i+1]<y[i]) {
	maxima.push_back(make_pair(x[i],y[i]));
      }
    }
69 70 71 72
    }
  }


73
  void transformMesh(Mesh *mesh, WorldVector<double> scale, WorldVector<double> shift, WorldVector<double> rotate)
Praetorius, Simon's avatar
Praetorius, Simon committed
74
  {
75

76 77 78
    int dow = Global::getGeo(WORLD);
    mtl::dense2D<double> Rx(dow,dow),Ry(dow,dow),Rz(dow,dow),T_(dow,dow),S(dow,dow);
    Rx = 1.0;
79 80
    Ry = 1.0;
    Rz = 1.0;
81 82 83 84
    S = 1.0;
    if (dow == 3) {
      Rx[1][1] = cos(rotate[0]);  Rx[1][2] = -sin(rotate[0]);
      Rx[2][1] = sin(rotate[0]);  Rx[2][2] = cos(rotate[0]);
85

86 87
      Ry[0][0] = cos(rotate[1]);  Ry[0][2] = sin(rotate[1]);
      Ry[2][0] = -sin(rotate[1]); Ry[2][2] = cos(rotate[1]);
88

89 90
      Rz[0][0] = cos(rotate[2]);  Rz[0][1] = -sin(rotate[2]);
      Rz[1][0] = sin(rotate[2]);  Rz[1][1] = cos(rotate[2]);
91

92 93 94 95
      S[0][0] = scale[0];  S[1][1] = scale[1];  S[2][2] = scale[2];
    } else if (dow == 2) {
      Rx[0][0] = cos(rotate[0]);  Rx[0][1] = -sin(rotate[0]);
      Rx[1][0] = sin(rotate[0]);  Rx[1][1] = cos(rotate[0]);
96

97 98
      Ry = 1.0;
      Rz = 1.0;
99

100 101 102
      S[0][0] = scale[0];  S[1][1] = scale[1];
    } else {
      S[0][0] = scale[0];
103 104
    }

105 106 107 108
    T_ = Rz*Ry*Rx*S;

    WorldMatrix<double> T;
    vector_operations::fillMatrix(T_,T);
109

110 111
    FixVec<WorldVector<double>, VERTEX> coords;
    deque<MacroElement*>::iterator macro;
112

113 114 115 116 117 118 119 120
    for(macro = mesh->firstMacroElement(); macro != mesh->endOfMacroElements(); macro++) {
      for (int i = 0; i < mesh->getDim()+1; ++i) {
	WorldVector<double> x = (*macro)->getCoord(i);
	WorldVector<double> x_temp = x;
	x.multMatrixVec(T,x_temp);
	x+= shift;
	(*macro)->setCoord(i, x);
      }
121 122 123
    }
  }

124 125
  
  void scaleMesh(Mesh *mesh, WorldVector<double> shift, WorldVector<double> scale)
Praetorius, Simon's avatar
Praetorius, Simon committed
126
  {
127 128 129 130 131 132 133 134 135 136 137
    deque<MacroElement*>::iterator macro;

    for(macro = mesh->firstMacroElement(); macro != mesh->endOfMacroElements(); macro++) {
      for (int i = 0; i < mesh->getDim()+1; ++i) {
	WorldVector<double> x = (*macro)->getCoord(i);
	x *= scale;
	x += shift;
	(*macro)->setCoord(i, x);
      }
    }
  }
138

139 140 141
  
  // scale by different values in all directions
  void scaleMesh(Mesh *mesh, WorldVector<double> scale)
Praetorius, Simon's avatar
Praetorius, Simon committed
142
  {
143 144 145 146
    WorldVector<double> shift; shift.set(0.0);
    scaleMesh(mesh,shift,scale);
  }
  
147

148 149
  // scale by different values in all directions
  void scaleMesh(std::vector<Mesh*> meshes, WorldVector<double> scale)
Praetorius, Simon's avatar
Praetorius, Simon committed
150
  {
151 152 153 154 155
    WorldVector<double> shift; shift.set(0.0);
    for (size_t i = 0; i < meshes.size(); i++)
      scaleMesh(meshes[i],shift,scale);
  }
  
156

157 158
  // scale and shift by the same values in all directions
  void scaleMesh(Mesh *mesh, double shift, double scale)
Praetorius, Simon's avatar
Praetorius, Simon committed
159
  {
160 161 162 163 164 165
    WorldVector<double> wShift, wScale;
    wShift.set(shift);
    wScale.set(scale);
    scaleMesh(mesh,wShift,wScale);
  }
  
166

167 168
  void read_dof_vector(const std::string file, DOFVector<double> *dofvec, long size)
  { FUNCNAME("Helpers::read_dof_vector()");
169

170 171 172
    ifstream in;
    string buf;
    long i;
173

174
    in.open(file.c_str());
Praetorius, Simon's avatar
Praetorius, Simon committed
175
    TEST_EXIT(in)("Could not read value file.\n");
176

177
    DOFIterator<double> dofIter(dofvec, USED_DOFS);
178

179 180 181 182 183 184 185
    getline(in, buf);
    while(buf.substr(0,14) != "vertex values:") getline(in, buf);
    for (i = 0, dofIter.reset(); !dofIter.end(); ++i, ++dofIter){
      in >> buf;
      (*dofIter) = atof(buf.c_str());
    }
    in.close();
186

Praetorius, Simon's avatar
Praetorius, Simon committed
187
    TEST_EXIT(i == size)("#vertices != dofIter.size()\n");
188
  }
189

190 191
  
  WorldVector<double> getMeshDimension(Mesh *mesh)
Praetorius, Simon's avatar
Praetorius, Simon committed
192
  {
193 194 195 196 197 198 199 200 201 202
    WorldVector<double> dimension; dimension.set(-1.e10);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int i = 0; i <= mesh->getDim(); i++) {
	WorldVector<double> &coords = elInfo->getMacroElement()->getCoord(i);
	for(int j = 0; j < coords.getSize(); ++j) {
	  dimension[j] = std::max(dimension[j], coords[j]);
	}
203
      }
204
      elInfo = stack.traverseNext(elInfo);
205
    }
206
    return dimension;
207
  }
208 209 210

  
  void getMeshDimension(Mesh *mesh, WorldVector<double> &min_corner, WorldVector<double> &max_corner)
Praetorius, Simon's avatar
Praetorius, Simon committed
211
  {
212 213 214 215 216 217 218 219 220 221 222 223
    max_corner.set(-1.e10);
    min_corner.set(1.e10);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL | Mesh::FILL_COORDS);
    while (elInfo) {
      for (int i = 0; i <= mesh->getDim(); i++) {
	WorldVector<double> &coords = elInfo->getMacroElement()->getCoord(i);
	for(int j = 0; j < coords.getSize(); ++j) {
	  max_corner[j] = std::max(max_corner[j], coords[j]);
	  min_corner[j] = std::min(min_corner[j], coords[j]);
	}
224
      }
225
      elInfo = stack.traverseNext(elInfo);
226 227 228 229
    }
  }

  
230 231 232 233 234
  void getNormalsWeighted(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals)
  {
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    int nBasisFcts = basisFcts->getNumber();
235
    
236
    std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
237
    
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
    DOFVector<double> numNormals(feSpace, "numNormals");
    numNormals.set(0.0);
      
    WorldVector<double> normal; normal.set(0.0);
    normals->set(normal);
    TraverseStack stack;
    ElInfo *elInfo = 
      stack.traverseFirst(mesh, -1, 
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
      
    while (elInfo) {
      elInfo->getElementNormal(normal);
      Element *el = elInfo->getElement();
      double det = elInfo->getDet();
      if (det <= DBL_TOL)
	det = elInfo->calcDet();
      
      basisFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices);
      
      for (int i = 0; i < nBasisFcts; i++) {
	(*normals)[localIndices[i]] += normal*(1.0/det);
      }
260

261 262 263 264 265 266 267 268
      elInfo = stack.traverseNext(elInfo);
    }
      
    DOFIterator<WorldVector<double> > normalsIter(normals, USED_DOFS); normalsIter.reset();
    for (; !normalsIter.end(); normalsIter++)
    {
      *normalsIter *= 1.0/norm(*normalsIter);
    }
269 270 271
  }

  
272 273 274 275 276 277 278 279
  void getNormals(FiniteElemSpace *feSpace, DOFVector<WorldVector<double> > *normals, DOFVector<WorldMatrix<double> > *gradNormals)
  {
    using namespace mtl;
    using namespace itl;
    
    Mesh *mesh = feSpace->getMesh();
    const BasisFunction *basisFcts = feSpace->getBasisFcts();
    size_t nBasisFcts = basisFcts->getNumber();
280

281 282
    std::vector<DegreeOfFreedom> localIndices(nBasisFcts);
    DOFVector<std::set<DegreeOfFreedom> > neighbors(feSpace, "neighbors");
283

284 285 286 287
    TraverseStack stack;
    ElInfo *elInfo =
      stack.traverseFirst(mesh, -1,
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
288

289 290
    while (elInfo) {
      basisFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
291

292 293 294 295
      for (size_t i = 0; i < nBasisFcts; i++) {
	for (size_t j = 0; j < nBasisFcts; j++) {
	  neighbors[localIndices[i]].insert(localIndices[j]);
	}
296
      }
297
      elInfo = stack.traverseNext(elInfo);
298 299
    }

300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384
    DOFVector<WorldVector<double> > coordsDOF(feSpace, "coords");
    mesh->getDofIndexCoords(coordsDOF);

    Helpers::getNormalsWeighted(feSpace, normals);
    VtkWriter::writeFile(normals, "normals_weighted.vtu");

  //   DOFVector<WorldMatrix<double> > *grdDOF = gradNormals;
  //   if (grdDOF == NULL)
  //     grdDOF = new DOFVector<WorldMatrix<double> >(feSpace, "grd");

    DOFIterator<WorldVector<double> > normalsIter(normals, USED_DOFS); normalsIter.reset();
    DOFIterator<WorldVector<double> > coordsIter(&coordsDOF, USED_DOFS); coordsIter.reset();
  /*  DOFIterator<WorldMatrix<double> > grdIter(grdDOF, USED_DOFS); grdIter.reset();*/
    DOFIterator<std::set<DegreeOfFreedom> > neighborsIter(&neighbors, USED_DOFS); neighborsIter.reset();
    for (; !neighborsIter.end(); neighborsIter++, normalsIter++, coordsIter++/*, grdIter++*/)
    {
      size_t numPoints = 10; // 6 oder 10
      std::set<DegreeOfFreedom> idx = *neighborsIter;
      std::set<DegreeOfFreedom>::iterator it, it2;

      for (size_t i = 1; i<3 && idx.size() < numPoints; i++) {
      for (it = idx.begin(); it != idx.end() && idx.size() < numPoints; it++) {
	for (it2 = (neighbors[*it]).begin(); it2 != (neighbors[*it]).end() && idx.size() < numPoints; it2++) {
	  idx.insert(*it2);
	}
      }
      }
      TEST_EXIT(idx.size()>=numPoints)("Nicht genügend Stützstellen gefunden!\n");
      size_t i = 0;
      mtl::dense2D<double> A(numPoints,6), B(6,6);
      mtl::dense_vector<double> b(numPoints, 1.0), c(6), y(6, 1.0);

      for (it = idx.begin(); it != idx.end() && i < numPoints; it++, i++) {
	WorldVector<double> x = coordsDOF[*it];
	A[i][0] = sqr(x[0]); A[i][1] = sqr(x[1]); A[i][2] = sqr(x[2]);
	A[i][3] = x[0]*x[1]; A[i][4] = x[0]*x[2]; A[i][5] = x[1]*x[2];
  //       if(numPoints>6)
  // 	A[i][6] = x[0]; A[i][7] = x[1]; A[i][8] = x[2]; A[i][9] = 1.0;
      }
      c = trans(A)*b;
      B = trans(A)*A;
  //     pc::ilu_0<mtl::dense2D<double> > P(B);
  //     basic_iteration<double>          iter(c, 100, 1.e-6);

  //     Solve Ax == b with left preconditioner P
  //     bicgstab(B, y, c, P, iter);

      y = mtl::matrix::lu_solve(B,c);

      WorldVector<double> oldNormal = (*normalsIter);

      (*normalsIter)[0] = 2.0*y[0]*(*coordsIter)[0] + y[3]*(*coordsIter)[1] + y[4]*(*coordsIter)[2]/* + (numPoints==10?y[6]:0.0)*/;
      (*normalsIter)[1] = 2.0*y[1]*(*coordsIter)[1] + y[3]*(*coordsIter)[0] + y[5]*(*coordsIter)[2]/* + (numPoints==10?y[7]:0.0)*/;
      (*normalsIter)[2] = 2.0*y[2]*(*coordsIter)[2] + y[4]*(*coordsIter)[0] + y[5]*(*coordsIter)[1]/* + (numPoints==10?y[8]:0.0)*/;

      double nrm = norm(*normalsIter);
      double diff = norm(oldNormal - *normalsIter);
      double signFactor = 1.0;
      if (diff > nrm)
	signFactor = -1.0;
      *normalsIter *= signFactor/nrm;

  //     if(gradNormals != NULL) {
  //       WorldMatrix<double> grad;
  //       grad[0][0] = 2.0*y[0]; grad[0][1] = y[3]; grad[0][2] = y[4];
  //       grad[1][0] = y[3]; grad[1][1] = 2.0*y[1]; grad[1][2] = y[5];
  //       grad[2][0] = y[4]; grad[2][1] = y[5]; grad[2][2] = 2.0*y[2];
  //       grad *= signFactor;
  //       
  //       double trace_ = 0.0;
  //       for (int i=0;i<3;i++)
  // 	trace_ += (grad[i][i]/nrm-sqr((*normalsIter)[i]));
  //       std::cout<<"nrm(n)="<<nrm<<std::endl;
  //       std::cout<<"trace(grad(n))="<<vector_operations::trace(grad)/nrm<<"="<<trace_<<std::endl;
  //       
  //       WorldMatrix<double> NN;
  //       NN.vecProduct(*normalsIter, *normalsIter);
  // 
  //       grad -= NN;
  //       grad *= 1.0/nrm;
  //       
  //       *grdIter = grad;
  // 
  //     }
    }
385

386 387 388
  //   if(gradNormals == NULL)
  //     delete grdDOF;
  }
389

390 391
  
  void getCurvature(DOFVector<WorldVector<double> >* normals,DOFVector<double>* curvature)
392
  {
393 394 395 396 397 398 399 400
    int dim = normals->getFeSpace()->getMesh()->getDim();
    int dow = Global::getGeo(WORLD);
    
    WorldVector<DOFVector<double>*> normalsVec;
    WorldVector<DOFVector<WorldVector<double> >*> gradNormals;
    for (int i = 0; i < dow; i++) {
      normalsVec[i] = new DOFVector<double>(normals->getFeSpace(), "n_i");
      gradNormals[i] = new DOFVector<WorldVector<double> >(normals->getFeSpace(), "grad_i");
401
    }
402 403 404 405 406 407 408 409 410 411 412 413
    transform(normals, &normalsVec);
    for (int i = 0; i < dow; i++)
      normalsVec[i]->getRecoveryGradient(gradNormals[i]);

    DOFIterator<WorldVector<double> > grd0Iter(gradNormals[0], USED_DOFS); grd0Iter.reset();
    DOFIterator<WorldVector<double> > grd1Iter(gradNormals[1], USED_DOFS); grd1Iter.reset();
    DOFIterator<WorldVector<double> > grd2Iter(gradNormals[2], USED_DOFS); grd2Iter.reset();
    DOFIterator<double> hIter(curvature, USED_DOFS); hIter.reset();
    
    for (; !grd0Iter.end(); grd0Iter++,grd1Iter++,grd2Iter++,hIter++)
    {
      *hIter = ((*grd0Iter)[0] + (*grd1Iter)[1] + (*grd2Iter)[2])/static_cast<double>(dim);
414
    }
415 416 417 418

    for (int i = 0; i < dow; i++) {
      delete normalsVec[i];
      delete gradNormals[i];
419 420 421 422
    }
  }

  
423 424
  // plots a vector with ascii-code
  void plot(std::vector<double> &values, int numRows, int numCols, std::string symbol)
425
  {
426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441
    unsigned num=values.size();
    unsigned cols=std::min(numCols,(int)num), rows=std::min(numRows,(int)num);
    unsigned step=(unsigned)floor(num/(double)cols);
    cols=(unsigned)ceil(num/(double)step);

    double maxVal=*max_element(values.begin(),values.end());
    double minVal=*min_element(values.begin(),values.end());
    std::stringstream ssMax; ssMax.setf(std::ios::fixed,std::ios::floatfield); ssMax.precision(2); ssMax<<maxVal;
    std::stringstream ssMin; ssMin.setf(std::ios::fixed,std::ios::floatfield); ssMin.precision(2); ssMin<<minVal;
    unsigned cols0 = std::max(ssMax.str().length(),ssMin.str().length());
    std::stringstream ssNum; ssNum<<num;

    std::string emptyLine(cols+cols0,' ');emptyLine.replace(cols0-1,1,"|");
    std::vector<std::string> lines(rows+1,emptyLine);

    double range=maxVal-minVal;
442
    
443 444 445 446 447 448 449
    for(unsigned i=0;i<cols; i+=1) {
      int v=(int)round((rows-1)*(values[i*step]-minVal)/range,0);
      int r=rows-v-1, c=i+cols0;
      
      lines[r].replace(c,1,symbol);
      lines[rows].replace(c,1,"-");
    }
450

451 452 453 454 455
    lines[0].replace(cols0-ssMax.str().length(),ssMax.str().length(),ssMax.str());
    lines[rows-1].replace(cols0-ssMin.str().length(),ssMin.str().length(),ssMin.str());
    lines[rows].replace(cols0-1,2," 0");
    lines[rows].erase(cols+cols0-1,1);
    lines[rows].append(ssNum.str());
456

457 458 459
    for(unsigned i=0;i<lines.size(); ++i)
      std::cout<<lines[i]<<std::endl;
  }
Praetorius, Simon's avatar
Helpers  
Praetorius, Simon committed
460

461
} // end namespace Helpers