Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64340528
Multiquadratics_2D_WENO.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, May 26, 05:47
Size
9 KB
Mime Type
text/x-c
Expires
Tue, May 28, 05:47 (2 d)
Engine
blob
Format
Raw Data
Handle
17886817
Attached To
R3721 ConSol
Multiquadratics_2D_WENO.cpp
View Options
//
// Created by Fabian Moenkeberg on 2019-12-16.
//
#include "Multiquadratics_2D_WENO.hpp"
#include "cblas.h"
Multiquadratics_2D_WENO
::
Multiquadratics_2D_WENO
(
int
order
,
double
eps
){
this
->
order
=
order
;
this
->
order_scheme
=
2
*
((
order
+
1
)
/
2
+
(
order
+
1
)
%
2
);
this
->
shapeParam
=
eps
;
this
->
dim
=
2
;
this
->
degPoly
=
order
-
1
;
this
->
name
=
MULTIQUADS
;
setMulti_Index
(
degPoly
+
1
);
}
Multiquadratics_2D_WENO
::
Multiquadratics_2D_WENO
(
int
order
,
int
degPoly
,
double
eps
){
this
->
order
=
1
;
// this->order_scheme = 2*((order+1)/2 + (order+1)%2);
this
->
order_scheme
=
order
+
1
;
this
->
shapeParam
=
eps
;
this
->
dim
=
2
;
this
->
degPoly
=
degPoly
;
this
->
name
=
MULTIQUADS
;
setMulti_Index
(
degPoly
+
1
);
}
std
::
vector
<
double
>
Multiquadratics_2D_WENO
::
evaluate
(
std
::
vector
<
std
::
vector
<
double
>>
evPoints
,
MeshBase
*
mesh
,
std
::
vector
<
int
>
stencil
,
std
::
vector
<
double
>
values
){
int
nStencils
=
(
int
)
stencil
.
size
();
std
::
vector
<
std
::
vector
<
std
::
vector
<
double
>>>
xB
(
nStencils
,
std
::
vector
<
std
::
vector
<
double
>>
(
NODE_NUM
,
std
::
vector
<
double
>
(
2
)));
std
::
vector
<
std
::
vector
<
double
>>
xi
(
NODE_NUM
,
std
::
vector
<
double
>
(
2
));
for
(
int
i
=
0
;
i
<
nStencils
;
i
++
){
xi
=
mesh
->
getXCoord
(
mesh
->
getCells
(
stencil
[
i
]));
xB
[
i
]
=
xi
;
}
double
rad
=
stableRad
/
mesh
->
getApproxStencilSize
(
stencil
);
// double rad = stableRad / sqrt(mesh->getSize(cellIndex));
std
::
vector
<
double
>
ret
=
Multiquadratics_2D_WENO
::
evaluateDirect
(
evPoints
,
xB
,
values
,
1
,
rad
);
return
ret
;
}
std
::
vector
<
double
>
Multiquadratics_2D_WENO
::
evaluateDirect
(
std
::
vector
<
std
::
vector
<
double
>>
evPoints
,
std
::
vector
<
std
::
vector
<
std
::
vector
<
double
>>>
xB
,
std
::
vector
<
double
>
values
,
double
eps
,
double
rad
){
double
eps0
=
1e-10
;
int
nEv
=
(
int
)
evPoints
.
size
();
std
::
vector
<
double
>
smInd
(
3
);
std
::
vector
<
std
::
vector
<
double
>>
ret
(
3
,
std
::
vector
<
double
>
(
nEv
,
0
));
std
::
vector
<
double
>
res
(
nEv
,
0
);
int
nRbf
=
(
int
)
xB
.
size
();
std
::
vector
<
double
>
coeff_old
;
int
nRbf0
=
std
::
min
((
int
)
xB
.
size
(),
3
);
int
nRbf_old
;
int
nP_old
;
std
::
vector
<
std
::
vector
<
std
::
vector
<
double
>>>
xB_
;
// std::vector<std::vector<std::vector<double>>> xB_.reserve(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
std
::
vector
<
double
>
values_
;
std
::
vector
<
double
>
values_old
;
int
degPolyMax_
=
std
::
max
(
getDegPoly
(
nRbf
),
0
);
for
(
int
degPoly_
=
0
;
degPoly_
<=
degPolyMax_
;
degPoly_
++
)
{
int
nP
=
multi_index
.
size
();
int
nP_tmp
=
0
;
int
sum_multi_index
=
multi_index
[
0
][
0
]
+
multi_index
[
0
][
1
];
for
(
int
i
=
0
;
i
<
nP
;
i
++
)
{
if
(
multi_index
[
i
][
0
]
+
multi_index
[
i
][
1
]
<=
degPoly_
)
{
nP_tmp
++
;
}
}
nP
=
nP_tmp
;
if
(
degPoly_
==
0
)
{
nRbf
=
1
;
xB_
.
resize
(
nRbf
,
std
::
vector
<
std
::
vector
<
double
>>
(
3
,
std
::
vector
<
double
>
(
2
)));
xB_
[
0
]
=
xB
[
0
];
values_
.
resize
(
nRbf
);
values_
[
0
]
=
values
[
0
];
}
else
if
(
degPoly_
==
1
)
{
nRbf
=
3
;
xB_
.
resize
(
nRbf
,
std
::
vector
<
std
::
vector
<
double
>>
(
3
,
std
::
vector
<
double
>
(
2
)));
values_
.
resize
(
nRbf
);
for
(
int
i
=
0
;
i
<
nRbf
;
i
++
){
xB_
[
i
]
=
xB
[
i
];
values_
[
i
]
=
values
[
i
];
}
}
else
{
nRbf
=
xB
.
size
();
// nRbf = 6;
xB_
.
resize
(
nRbf
,
std
::
vector
<
std
::
vector
<
double
>>
(
3
,
std
::
vector
<
double
>
(
2
)));
values_
.
resize
(
nRbf
);
for
(
int
i
=
0
;
i
<
nRbf
;
i
++
){
xB_
[
i
]
=
xB
[
i
];
values_
[
i
]
=
values
[
i
];
}
}
int
m
=
nRbf
+
nP
;
eps
=
eps
*
rad
;
int
approxOrder
=
std
::
min
(
degPoly_
+
1
,
order
);
getInterpCoeff
(
xB_
,
values_
,
eps
,
nP
,
approxOrder
);
std
::
vector
<
double
>
values_tmp
=
values_
;
if
(
nRbf
>
1
){
for
(
int
i
=
0
;
i
<
nRbf_old
;
i
++
){
values_tmp
[
i
]
-=
values_old
[
i
];
}
for
(
int
i
=
0
;
i
<
nP_old
;
i
++
){
values_tmp
[
i
+
nRbf
]
-=
values_old
[
i
+
nRbf_old
];
}
}
values_old
=
values_
;
std
::
vector
<
std
::
vector
<
double
>>
EvA0
(
nEv
,
std
::
vector
<
double
>
(
nRbf
));
MVEvaluationMatrix
(
evPoints
,
xB_
,
eps
,
EvA0
,
approxOrder
);
std
::
vector
<
std
::
vector
<
double
>>
PA
(
nEv
,
std
::
vector
<
double
>
(
nP
));
MVPolynomialEvaluationMatrixV
(
evPoints
,
xB_
,
PA
,
eps
,
nP
);
double
*
EvA
=
new
double
[
nEv
*
m
];
for
(
int
i
=
0
;
i
<
nEv
;
i
++
)
{
for
(
int
j
=
0
;
j
<
nRbf
;
j
++
)
{
EvA
[
i
+
j
*
nEv
]
=
EvA0
[
i
][
j
];
}
for
(
int
j
=
0
;
j
<
nP
;
j
++
)
{
EvA
[
i
+
(
j
+
nRbf
)
*
nEv
]
=
PA
[
i
][
j
];
}
}
if
(
nRbf
>
1
){
smInd
[
degPoly_
]
=
1.0
/
(
ptwiseSmInd
(
xB
[
0
],
xB_
,
values_tmp
,
eps
)
+
eps0
)
+
1
;
}
else
{
// double a1,a2,a3,b1,b2,b3;
// if (degPolyMax_>0){
// for (int i = 0; i < 3; i++){
//
// }
// }
smInd
[
0
]
=
1.0
/
eps0
+
1
;
// double tmp = 0;
// for (int jj = 1; jj < nRbf0; jj++){
// tmp += (values[jj] - values[0])*(values[jj] - values[0]);
// }
// smInd[0] = 1.0/(tmp*std::sqrt(calcSize(xB[0])) + eps0)+1;
}
double
*
sVecPt
=
values_
.
data
();
double
*
y
=
new
double
[
nEv
];
cblas_dgemv
(
CblasColMajor
,
CBLAS_TRANSPOSE
::
CblasNoTrans
,
(
int
)
nEv
,
(
int
)
(
nRbf
+
nP
),
1.0
,
EvA
,
(
int
)
nEv
,
sVecPt
,
1
,
0.0
,
y
,
(
int
)
1
);
for
(
int
i
=
0
;
i
<
nEv
;
i
++
){
ret
[
degPoly_
][
i
]
=
y
[
i
];
}
// if (nRbf == 6){
// getInterpCoeff(xB_, values, eps, nP, approxOrder);
// double *sVecPt = values.data();
// cblas_dgemv(CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, (int) nEv, (int) (nRbf + nP), 1.0, EvA, (int) nEv,
// sVecPt, 1, 0.0, y, (int) 1);
//
// for (int i = 0; i < nEv; i++){
// for ( int j = 1; j < degPolyMax_; j++) {
// ret[0][i] += ret[j][i];
// }
// ret[1][i] = y[i];
// }
// }
// values_old = values_;
nP_old
=
nP
;
nRbf_old
=
nRbf
;
delete
[]
y
;
delete
[]
EvA
;
// delete[] sVecPt;
}
double
sum_smInd
=
0
;
for
(
int
i
=
0
;
i
<=
degPolyMax_
;
i
++
){
sum_smInd
+=
smInd
[
i
];
}
for
(
int
j
=
0
;
j
<
degPolyMax_
;
j
++
){
for
(
int
i
=
0
;
i
<
nEv
;
i
++
){
res
[
i
]
+=
(
ret
[
degPolyMax_
-
j
][
i
]
-
ret
[
degPolyMax_
-
j
-
1
][
i
])
*
smInd
[
degPolyMax_
-
j
]
/
sum_smInd
*
(
degPolyMax_
+
1
);
// ret[degPolyMax_-j][i] += ret[degPolyMax_-j][i]*smInd[degPolyMax_-j]/sum_smInd;
}
}
for
(
int
i
=
0
;
i
<
nEv
;
i
++
){
res
[
i
]
+=
ret
[
0
][
i
]
*
smInd
[
0
]
/
sum_smInd
*
(
degPolyMax_
+
1
);
}
std
::
cout
<<
smInd
[
0
]
/
sum_smInd
<<
", "
<<
smInd
[
1
]
/
sum_smInd
<<
", "
<<
smInd
[
2
]
/
sum_smInd
<<
"
\n
"
;
return
res
;
}
double
Multiquadratics_2D_WENO
::
ptwiseSmInd
(
std
::
vector
<
std
::
vector
<
double
>>
xBi
,
std
::
vector
<
std
::
vector
<
std
::
vector
<
double
>>>
xB
,
std
::
vector
<
double
>
values
,
double
eps
){
int
nRbf
=
(
int
)
xB
.
size
();
int
degPoly_
=
std
::
max
(
getDegPoly
(
nRbf
),
0
);
double
Ind
=
0
;
// double s = 0;
// double p = 1;
// double dist = 0;
std
::
vector
<
double
>
xM
(
2
,
0.0
);
xM
[
0
]
=
(
xBi
[
0
][
0
]
+
xBi
[
1
][
0
]
+
xBi
[
2
][
0
])
/
NODE_NUM
;
xM
[
1
]
=
(
xBi
[
0
][
1
]
+
xBi
[
1
][
1
]
+
xBi
[
2
][
1
])
/
NODE_NUM
;
// double delta;
// double d_0, d_1;
for
(
int
i
=
0
;
i
<
nRbf
;
i
++
){
// Ind += std::pow(values[i],2);
// delta = eps*std::sqrt(calcSize(xB[i]));
Ind
+=
values
[
i
]
*
values
[
i
];
// p *= delta*delta;
// s += delta*delta;
// d_0 = xM[0]-calcMidPt(xB[i],0);
// d_1 = xM[1]-calcMidPt(xB[i],1);
// dist += eps*std::sqrt(d_0*d_0 + d_1*d_1);
}
// Ind = Ind*p*p/(s*s);
Ind
=
0
;
Ind
/=
nRbf
;
double
delta0
=
std
::
sqrt
(
calcSize
(
xBi
));
for
(
int
j
=
1
;
j
<=
degPoly_
;
j
++
){
// double Ds = 0;
// for (int i = 0; i < nRbf;i++){
// Ds += values[i]*DF(std::vector<double>{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j);
// }
// Ds += DPolynomial(xM, xB, values, j, values.size() - nRbf, eps);
// Ind += std::pow(Ds,2);
for
(
int
j2
=
0
;
j2
<=
j
;
j2
++
){
double
Ds
=
0
;
for
(
int
i
=
0
;
i
<
nRbf
;
i
++
){
Ds
+=
values
[
i
]
*
DF
(
std
::
vector
<
double
>
{
eps
*
(
xM
[
0
]
-
calcMidPt
(
xB
[
i
],
0
)),
eps
*
(
xM
[
1
]
-
calcMidPt
(
xB
[
i
],
1
))},
j2
,
j
-
j2
,
eps
);
}
// Ds += DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
Ind
+=
std
::
pow
(
Ds
,
2
)
*
my_pow
(
delta0
,
j
);
}
}
// Ind *= dist;
return
std
::
fabs
(
Ind
);
}
int
Multiquadratics_2D_WENO
::
getDegPoly
(
int
n
){
int
degPoly_
;
if
(
false
){
degPoly_
=
1
;
}
else
{
degPoly_
=
(
int
)
floor
(
-
1.5
+
sqrt
(
1
+
8
*
(
n
))
/
2
);
}
return
degPoly_
;
}
Event Timeline
Log In to Comment