Peano 4
Loading...
Searching...
No Matches
TP_PunctureTracker.h
Go to the documentation of this file.
1#include <iostream>
2#include <fstream>
3#include <limits>
4#include <string>
5#include <cstring>
6
7#include "tarch/la/Vector.h"
8
9#pragma once
10
11
12inline std::istream& ignoreline(std::fstream& in, std::ifstream::pos_type& pos)
13{
14 pos = in.tellg();
15 return in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
16}
17
18inline std::string getLastLine(std::fstream& in)
19{
20 std::ifstream::pos_type pos = in.tellg();
21
22 std::ifstream::pos_type lastPos;
23 while (in >> std::ws && ignoreline(in, lastPos))
24 pos = lastPos;
25
26 in.clear();
27 in.seekg(pos);
28
29 std::string line;
30 std::getline(in, line);
31 return line;
32}
33
34inline double linearInter(double x1, double f1, double x2, double f2, double target){
35 if ((x1-x2)<10e-8){
36 return (f1+f2)/2;
37 } else {
38 return f2*((target-x1)/(x2-x1))+f1*((x2-target)/(x2-x1));
39 }
40}
41
42
43inline void CoorReadIn(double* coor, std::string line)
44{
45 int index=0;
46 std::string var[5]={"","","","",""};
47 for (auto x : line){
48 if (x == ' '){
49 coor[index]=std::stod(var[index]);
50 index++;
51 } else {
52 var[index]=var[index]+x;
53 }
54 }
55 coor[index]=std::stod(var[index]);
56}
57
58inline void ConsReadIn(double* cons, std::string line)
59{
60 int index=0;
61 std::string var[7]={"","","","","","",""};
62 for (auto x : line){
63 if (x == ' '){
64 cons[index]=std::stod(var[index]);
65 index++;
66 } else {
67 var[index]=var[index]+x;
68 }
69 }
70 cons[index]=std::stod(var[index]);
71}
72
73
75 const double* coor,
77 double volumeH,
78 int patchSize
79) {
80 //std::cout<<Offset(0)<<" "<<Offset(1)<<" "<<Offset(2)<<" "<< volumeH <<std::endl;
81 tarch::la::Vector<Dimensions*2,int> index= {0,0,0,0,0,0};
82 for (int i=0;i<patchSize;i++){
83 if (coor[0]>(Offset(0)+i*volumeH) and coor[0]<(Offset(0)+(i+1)*volumeH))
84 {index(0)=i;
85 if (coor[0]<(Offset(0)+(i+0.5)*volumeH))
86 {index(3)=-1;}
87 else if (coor[0]>(Offset(0)+(i+0.5)*volumeH))
88 {index(3)=1;}
89 }
90 }
91 index(0)+=3; //to include the extra layer
92 //std::cout<<coor[0]<<" "<<(Offset(0)+index(0)*volumeH)<<" "<<(Offset(0)+(index(0)+1)*volumeH)<<" "<< (Offset(0)+(index(0)+0.5)*volumeH) <<std::endl;
93 for (int i=0;i<patchSize;i++){
94 if (coor[1]>(Offset(1)+i*volumeH) and coor[1]<(Offset(1)+(i+1)*volumeH))
95 {index(1)=i;
96 if (coor[1]<(Offset(1)+(i+0.5)*volumeH))
97 {index(4)=-1;}
98 else if (coor[1]>(Offset(1)+(i+0.5)*volumeH))
99 {index(4)=1;}
100 }
101 }
102 index(1)+=3; //to include the extra layer
103 //std::cout<<coor[1]<<" "<<(Offset(1)+index(1)*volumeH)<<" "<<(Offset(1)+(index(1)+1)*volumeH)<<" "<< (Offset(1)+(index(1)+0.5)*volumeH) <<std::endl;
104 for (int i=0;i<patchSize;i++){
105 if (coor[2]>(Offset(2)+i*volumeH) and coor[2]<(Offset(2)+(i+1)*volumeH))
106 {index(2)=i;
107 if (coor[2]<(Offset(2)+(i+0.5)*volumeH))
108 {index(5)=-1;}
109 else if (coor[2]>(Offset(2)+(i+0.5)*volumeH))
110 {index(5)=1;}
111 }
112 }
113 index(2)+=3; //to include the extra layer
114 return index;
115}
116
117//output index for 8 nearest cells
119{
120 int x[2],y[2],z[2];
121
122 if (IndexOfCell(3)==1) {
123 x[0]=IndexOfCell(0);x[1]=IndexOfCell(0)+1;
124 } else if (IndexOfCell(3)==-1) {
125 x[0]=IndexOfCell(0)-1;x[1]=IndexOfCell(0);
126 } else {
127 x[0]=IndexOfCell(0);x[1]=IndexOfCell(0);}
128
129 if (IndexOfCell(4)==1) {
130 y[0]=IndexOfCell(1);y[1]=IndexOfCell(1)+1;
131 } else if (IndexOfCell(4)==-1) {
132 y[0]=IndexOfCell(1)-1;y[1]=IndexOfCell(1);
133 } else {
134 y[0]=IndexOfCell(1);y[1]=IndexOfCell(1);}
135
136 if (IndexOfCell(5)==1) {
137 z[0]=IndexOfCell(2);z[1]=IndexOfCell(2)+1;
138 } else if (IndexOfCell(5)==-1) {
139 z[0]=IndexOfCell(2)-1;z[1]=IndexOfCell(2);
140 } else {
141 z[0]=IndexOfCell(2);z[1]=IndexOfCell(2);}
142
143 for(int i=0;i<2;i++)
144 for(int j=0;j<2;j++)
145 for(int k=0;k<2;k++){
146 InterIndex[i*4+j*2+k]={x[i],y[j],z[k]};
147 }
148
149}
150
151
152//do the real interpolation
153inline void Interpolation(
154 double* result,
156 double* raw,
157 const double* coor,
159 double volumeH,
160 int patchSize
161){
162 int inter_number=4;
163 //calculate the actual coordinates
164 double CoorsForInter1[8][inter_number];
165 double raw1[8][inter_number];
166 for(int i=0;i<2;i++)
167 for(int j=0;j<2;j++)
168 for(int k=0;k<2;k++){
169 for (int m=0;m<inter_number;m++){
170 CoorsForInter1[i*4+j*2+k][m]=Offset(m)+(IndexForInter[i*4+j*2+k](m)-0.5)*volumeH;
171 raw1[i*4+j*2+k][m]=raw[(i*4+j*2+k)*inter_number+m];
172 }
173 }
174
175 //first interpolate along x axis
176 double CoorsForInter2[4][inter_number];
177 double raw2[4][inter_number];
178 for (int n=0;n<4;n++){
179 CoorsForInter2[n][0]=coor[0];
180 CoorsForInter2[n][1]=CoorsForInter1[n][1];
181 CoorsForInter2[n][2]=CoorsForInter1[n][2];
182 for (int m=0;m<inter_number;m++){
183 raw2[n][m]=linearInter(CoorsForInter1[n][0],raw1[n][m],CoorsForInter1[n+4][0],raw1[n+4][m],coor[0]);
184 }
185 }
186
187 //second interpolate along y axis
188 double CoorsForInter3[2][inter_number];
189 double raw3[2][inter_number];
190 for (int n=0;n<2;n++){
191 CoorsForInter3[n][0]=coor[0];
192 CoorsForInter3[n][1]=coor[1];
193 CoorsForInter3[n][2]=CoorsForInter1[n][2];
194 for (int m=0;m<inter_number;m++){
195 raw3[n][m]=linearInter(CoorsForInter2[n][1],raw2[n][m],CoorsForInter2[n+2][1],raw2[n+2][m],coor[1]);
196 }
197 }
198
199 //finally interpolate along z axis
200 double CoorsForInter4[1][inter_number];
201 double raw4[1][inter_number];
202 for (int n=0;n<1;n++){
203 CoorsForInter4[n][0]=coor[0];
204 CoorsForInter4[n][1]=coor[1];
205 CoorsForInter4[n][2]=coor[2];
206 for (int m=0;m<inter_number;m++){
207 raw4[n][m]=linearInter(CoorsForInter3[n][2],raw3[n][m],CoorsForInter3[n+1][2],raw3[n+1][m],coor[2]);
208 }
209 }
210
211 result[0]=raw4[0][0]; result[1]=raw4[0][1]; result[2]=raw4[0][2];result[3]=raw4[0][3];
212}
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
tarch::la::Vector< Dimensions *2, int > FindCellIndex(const double *coor, tarch::la::Vector< Dimensions, double > Offset, double volumeH, int patchSize)
std::string getLastLine(std::fstream &in)
void CoorReadIn(double *coor, std::string line)
void FindInterIndex(tarch::la::Vector< Dimensions, int > *InterIndex, tarch::la::Vector< Dimensions *2, int > IndexOfCell)
void Interpolation(double *result, tarch::la::Vector< Dimensions, int > *IndexForInter, double *raw, const double *coor, tarch::la::Vector< Dimensions, double > Offset, double volumeH, int patchSize)
double linearInter(double x1, double f1, double x2, double f2, double target)
void ConsReadIn(double *cons, std::string line)
std::istream & ignoreline(std::fstream &in, std::ifstream::pos_type &pos)
Simple vector class.
Definition Vector.h:134