|
|
@@ -53,6 +53,61 @@ template <class Real, Integer DIM, Integer ORDER> class Basis {
|
|
|
}
|
|
|
return nodes_;
|
|
|
}
|
|
|
+ static const Vector<ValueType>& QuadWts() {
|
|
|
+ static Vector<ValueType> wts(Size());
|
|
|
+ { // Set nodes_
|
|
|
+ static std::mutex mutex;
|
|
|
+ static std::atomic<Integer> first_time(true);
|
|
|
+ if (first_time.load(std::memory_order_relaxed)) {
|
|
|
+ std::lock_guard<std::mutex> guard(mutex);
|
|
|
+ if (first_time.load(std::memory_order_relaxed)) {
|
|
|
+ StaticArray<ValueType,ORDER> wts_1d;
|
|
|
+ { // Set wts_1d
|
|
|
+ Vector<ValueType> x_(ORDER);
|
|
|
+ ChebBasis<ValueType>::template Nodes<1>(ORDER, x_);
|
|
|
+
|
|
|
+ Vector<ValueType> V_cheb(ORDER * ORDER);
|
|
|
+ { // Set V_cheb
|
|
|
+ Vector<ValueType> I(ORDER*ORDER);
|
|
|
+ I = 0;
|
|
|
+ for (Long i = 0; i < ORDER; i++) I[i*ORDER+i] = 1;
|
|
|
+ ChebBasis<ValueType>::template Approx<1>(ORDER, I, V_cheb);
|
|
|
+ }
|
|
|
+ Matrix<ValueType> M(ORDER, ORDER, V_cheb.begin());
|
|
|
+
|
|
|
+ Vector<ValueType> w_sample(ORDER);
|
|
|
+ for (Integer i = 0; i < ORDER; i++) {
|
|
|
+ w_sample[i] = (i % 2 ? 0 : -(ORDER/(ValueType)(i*i-1)));
|
|
|
+ }
|
|
|
+ for (Integer j = 0; j < ORDER; j++) {
|
|
|
+ wts_1d[j] = 0;
|
|
|
+ for (Integer i = 0; i < ORDER; i++) {
|
|
|
+ wts_1d[j] += M[j][i] * w_sample[i] / ORDER;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ wts[0] = 1;
|
|
|
+ Integer N = 1;
|
|
|
+ for (Integer d = 0; d < DIM; d++) {
|
|
|
+ for (Integer j = 1; j < ORDER; j++) {
|
|
|
+ for (Integer i = 0; i < N; i++) {
|
|
|
+ wts[j*N+i] = wts[i] * wts_1d[j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ for (Integer i = 0; i < N; i++) {
|
|
|
+ wts[i] *= wts_1d[0];
|
|
|
+ }
|
|
|
+ N *= ORDER;
|
|
|
+ }
|
|
|
+
|
|
|
+ std::atomic_thread_fence(std::memory_order_seq_cst);
|
|
|
+ first_time.store(false);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return wts;
|
|
|
+ }
|
|
|
|
|
|
static void Grad(Vector<Basis>& dX, const Vector<Basis>& X) {
|
|
|
static Matrix<ValueType> GradOp[DIM];
|
|
|
@@ -119,7 +174,7 @@ template <class Real, Integer DIM, Integer ORDER> class Basis {
|
|
|
Matrix<ValueType> M(Size(), N);
|
|
|
{ // Set M
|
|
|
auto nodes = Basis<ValueType,1,ORDER>::Nodes();
|
|
|
- Integer NN = nodes.Dim(1);
|
|
|
+ Integer NN = Basis<ValueType,1,ORDER>::Size();
|
|
|
Matrix<ValueType> M_(NN, DIM*N);
|
|
|
for (Long i = 0; i < DIM*N; i++) {
|
|
|
ValueType x = X[0][i];
|
|
|
@@ -171,6 +226,45 @@ template <class Real, Integer DIM, Integer ORDER> class Basis {
|
|
|
}
|
|
|
}
|
|
|
|
|
|
+ Basis operator+(Basis X) const {
|
|
|
+ for (Long i = 0; i < Size(); i++) X[i] = (*this)[i] + X[i];
|
|
|
+ return X;
|
|
|
+ }
|
|
|
+ Basis operator-(Basis X) const {
|
|
|
+ for (Long i = 0; i < Size(); i++) X[i] = (*this)[i] - X[i];
|
|
|
+ return X;
|
|
|
+ }
|
|
|
+ Basis operator*(Basis X) const {
|
|
|
+ for (Long i = 0; i < Size(); i++) X[i] = (*this)[i] * X[i];
|
|
|
+ return X;
|
|
|
+ }
|
|
|
+ Basis operator*(Real a) const {
|
|
|
+ Basis X = (*this);
|
|
|
+ for (Long i = 0; i < Size(); i++) X[i] *= a;
|
|
|
+ return X;
|
|
|
+ }
|
|
|
+ Basis& operator+=(const Basis& X) {
|
|
|
+ for (Long i = 0; i < Size(); i++) (*this)[i] += X[i];
|
|
|
+ return *this;
|
|
|
+ }
|
|
|
+ Basis& operator-=(const Basis& X) {
|
|
|
+ for (Long i = 0; i < Size(); i++) (*this)[i] -= X[i];
|
|
|
+ return *this;
|
|
|
+ }
|
|
|
+ Basis& operator*=(const Basis& X) {
|
|
|
+ for (Long i = 0; i < Size(); i++) (*this)[i] *= X[i];
|
|
|
+ return *this;
|
|
|
+ }
|
|
|
+ Basis& operator*=(Real a) {
|
|
|
+ for (Long i = 0; i < Size(); i++) (*this)[i] *= a;
|
|
|
+ return *this;
|
|
|
+ }
|
|
|
+ Basis& operator=(Real a) {
|
|
|
+ for (Long i = 0; i < Size(); i++) (*this)[i] = a;
|
|
|
+ return *this;
|
|
|
+ }
|
|
|
+
|
|
|
+
|
|
|
const ValueType& operator[](Long i) const {
|
|
|
SCTL_ASSERT(i < Size());
|
|
|
return NodeValues_[i];
|
|
|
@@ -368,7 +462,7 @@ template <class Real> class Quadrature {
|
|
|
|
|
|
|
|
|
|
|
|
- template <class DensityBasis, class ElemList, class Kernel> static void SetupSingular(Matrix<Real>& M_singular, const Matrix<Real>& trg_nds, const ElemList& elem_lst, const Kernel& kernel, Integer order_singular = 10, Integer order_direct = 10) {
|
|
|
+ template <class DensityBasis, class ElemList, class Kernel> static void SetupSingular(Matrix<Real>& M_singular, const Matrix<Real>& trg_nds, const ElemList& elem_lst, const Kernel& kernel, Integer order_singular = 10, Integer order_direct = 10, Real Rqbx = 0) {
|
|
|
using CoordBasis = typename ElemList::CoordBasis;
|
|
|
using CoordEvalOpType = typename CoordBasis::EvalOpType;
|
|
|
using DensityEvalOpType = typename DensityBasis::EvalOpType;
|
|
|
@@ -382,17 +476,43 @@ template <class Real> class Quadrature {
|
|
|
const Integer Ntrg = trg_nds.Dim(1);
|
|
|
SCTL_ASSERT(trg_nds.Dim(0) == ElemDim);
|
|
|
|
|
|
- Vector<Real> Xt;
|
|
|
- { // Set Xt
|
|
|
- auto Meval = CoordBasis::SetupEval(trg_nds);
|
|
|
- eval_basis(Xt, elem_lst.ElemVector(), ElemList::CoordDim(), trg_nds.Dim(1), Meval);
|
|
|
- }
|
|
|
- SCTL_ASSERT(Xt.Dim() == Nelem * Ntrg * CoordDim);
|
|
|
-
|
|
|
const Vector<CoordBasis>& X = elem_lst.ElemVector();
|
|
|
Vector<CoordBasis> dX;
|
|
|
CoordBasis::Grad(dX, X);
|
|
|
|
|
|
+ Vector<Real> Xt, Xnt;
|
|
|
+ { // Set Xt, Xnt
|
|
|
+ auto Meval = CoordBasis::SetupEval(trg_nds);
|
|
|
+ eval_basis(Xt, X, CoordDim, trg_nds.Dim(1), Meval);
|
|
|
+
|
|
|
+ Xnt = Xt;
|
|
|
+ Vector<Real> dX_;
|
|
|
+ eval_basis(dX_, dX, 2*CoordDim, trg_nds.Dim(1), Meval);
|
|
|
+
|
|
|
+ for (Long i = 0; i < Ntrg; i++) {
|
|
|
+ for (Long j = 0; j < Nelem; j++) {
|
|
|
+ auto Xn = Xnt.begin() + (j*Ntrg+i)*CoordDim;
|
|
|
+ auto dX0 = dX_.begin() + (j*Ntrg+i)*2*CoordDim;
|
|
|
+
|
|
|
+ StaticArray<Real,CoordDim> normal;
|
|
|
+ normal[0] = dX0[2]*dX0[5] - dX0[4]*dX0[3];
|
|
|
+ normal[1] = dX0[4]*dX0[1] - dX0[0]*dX0[5];
|
|
|
+ normal[2] = dX0[0]*dX0[3] - dX0[2]*dX0[1];
|
|
|
+ Real Xa = sctl::sqrt<Real>(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
|
|
|
+ Real invXa = 1/Xa;
|
|
|
+ normal[0] *= invXa;
|
|
|
+ normal[1] *= invXa;
|
|
|
+ normal[2] *= invXa;
|
|
|
+
|
|
|
+ Real sqrt_Xa = sqrt<Real>(Xa);
|
|
|
+ Xn[0] = normal[0]*sqrt_Xa*Rqbx;
|
|
|
+ Xn[1] = normal[1]*sqrt_Xa*Rqbx;
|
|
|
+ Xn[2] = normal[2]*sqrt_Xa*Rqbx;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ SCTL_ASSERT(Xt.Dim() == Nelem * Ntrg * CoordDim);
|
|
|
+
|
|
|
auto& M = M_singular;
|
|
|
M.ReInit(Nelem * KDIM0 * DensityBasis::Size(), KDIM1 * Ntrg);
|
|
|
#pragma omp parallel for schedule(static)
|
|
|
@@ -405,7 +525,7 @@ template <class Real> class Quadrature {
|
|
|
trg_node_[k] = trg_nds[k][i];
|
|
|
}
|
|
|
Vector<Real> trg_node(ElemDim, trg_node_, false);
|
|
|
- DuffyQuad<ElemDim>(quad_nds, quad_wts, trg_node, order_singular);
|
|
|
+ DuffyQuad<ElemDim>(quad_nds, quad_wts, trg_node, order_singular, fabs(Rqbx));
|
|
|
}
|
|
|
const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds);
|
|
|
Integer Nnds = quad_wts.Dim();
|
|
|
@@ -441,12 +561,33 @@ template <class Real> class Quadrature {
|
|
|
|
|
|
for (Long j = 0; j < Nelem; j++) {
|
|
|
Matrix<Real> M__(Nnds * KDIM0, KDIM1);
|
|
|
- { // Set kernel matrix M__
|
|
|
- const Vector<Real> X0_(CoordDim, (Iterator<Real>)Xt.begin() + (j * Ntrg + i) * CoordDim, false);
|
|
|
+ if (Rqbx == 0) { // Set kernel matrix M__
|
|
|
+ const Vector<Real> X0_(CoordDim, Xt.begin() + (j * Ntrg + i) * CoordDim, false);
|
|
|
const Vector<Real> X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false);
|
|
|
const Vector<Real> Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false);
|
|
|
kernel.template KernelMatrix<Real>(M__, X0_, X__, Xn__);
|
|
|
+ } else {
|
|
|
+ Vector<Real> X0_(CoordDim);
|
|
|
+ constexpr Integer qbx_order = 6;
|
|
|
+ StaticArray<Matrix<Real>,qbx_order> M___;
|
|
|
+ for (Integer k = 0; k < qbx_order; k++) { // Set kernel matrix M___
|
|
|
+ for (Integer kk = 0; kk < CoordDim; kk++) X0_[kk] = Xt[(j * Ntrg + i) * CoordDim + kk] + (k+1) * Xnt[(j * Ntrg + i) * CoordDim + kk];
|
|
|
+ const Vector<Real> X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false);
|
|
|
+ const Vector<Real> Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false);
|
|
|
+ kernel.template KernelMatrix<Real>(M___[k], X0_, X__, Xn__);
|
|
|
+ }
|
|
|
+
|
|
|
+ for (Long k = 0; k < Nnds * KDIM0 * KDIM1; k++) {
|
|
|
+ M__[0][k] = 0;
|
|
|
+ M__[0][k] += 6*M___[0][0][k];
|
|
|
+ M__[0][k] += -15*M___[1][0][k];
|
|
|
+ M__[0][k] += 20*M___[2][0][k];
|
|
|
+ M__[0][k] += -15*M___[3][0][k];
|
|
|
+ M__[0][k] += 6*M___[4][0][k];
|
|
|
+ M__[0][k] += -1*M___[5][0][k];
|
|
|
+ }
|
|
|
}
|
|
|
+
|
|
|
for (Long k0 = 0; k0 < KDIM0; k0++) {
|
|
|
for (Long k1 = 0; k1 < KDIM1; k1++) {
|
|
|
for (Long l = 0; l < DensityBasis::Size(); l++) {
|
|
|
@@ -1089,8 +1230,8 @@ template <class Real> class Quadrature {
|
|
|
for (Integer i = 0; i < 2; i++) { // iterate
|
|
|
Matrix<Real> X_, dX_;
|
|
|
for (Integer k = 0; k < ElemDim; k++) {
|
|
|
- u0(k,0) = std::min(1.0, u0(k,0));
|
|
|
- u0(k,0) = std::max(0.0, u0(k,0));
|
|
|
+ u0(k,0) = std::min<Real>(1.0, u0(k,0));
|
|
|
+ u0(k,0) = std::max<Real>(0.0, u0(k,0));
|
|
|
}
|
|
|
const auto eval_op = CoordBasis::SetupEval(Matrix<Real>(ElemDim,1,u0.begin(),false));
|
|
|
CoordBasis::Eval(X_, Vector<CoordBasis>(CoordDim,(Iterator<CoordBasis>)X.begin()+src_idx*CoordDim,false),eval_op);
|
|
|
@@ -1419,6 +1560,11 @@ template <class Real> class Quadrature {
|
|
|
public:
|
|
|
|
|
|
template <class DensityBasis, class ElemList, class Kernel> void Setup(const ElemList& elem_lst, const Vector<Real>& Xt, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) {
|
|
|
+ Xt_.ReInit(0);
|
|
|
+ M_singular.ReInit(0,0);
|
|
|
+ M_near_singular.ReInit(0,0);
|
|
|
+ pair_lst.ReInit(0);
|
|
|
+
|
|
|
order_direct_ = order_direct;
|
|
|
period_length_ = period_length;
|
|
|
comm_ = comm;
|
|
|
@@ -1438,7 +1584,12 @@ template <class Real> class Quadrature {
|
|
|
Profile::Toc();
|
|
|
}
|
|
|
|
|
|
- template <class DensityBasis, class PotentialBasis, class ElemList, class Kernel> void Setup(const ElemList& elem_lst, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) {
|
|
|
+ template <class DensityBasis, class PotentialBasis, class ElemList, class Kernel> void Setup(const ElemList& elem_lst, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm, Real Rqbx = 0) {
|
|
|
+ Xt_.ReInit(0);
|
|
|
+ M_singular.ReInit(0,0);
|
|
|
+ M_near_singular.ReInit(0,0);
|
|
|
+ pair_lst.ReInit(0);
|
|
|
+
|
|
|
order_direct_ = order_direct;
|
|
|
period_length_ = period_length;
|
|
|
comm_ = comm;
|
|
|
@@ -1475,7 +1626,7 @@ template <class Real> class Quadrature {
|
|
|
}
|
|
|
|
|
|
Profile::Tic("SetupSingular", &comm_);
|
|
|
- SetupSingular<DensityBasis>(M_singular, PotentialBasis::Nodes(), elem_lst, kernel, order_singular, order_direct_);
|
|
|
+ SetupSingular<DensityBasis>(M_singular, PotentialBasis::Nodes(), elem_lst, kernel, order_singular, order_direct_, Rqbx);
|
|
|
Profile::Toc();
|
|
|
|
|
|
Profile::Tic("SetupNearSingular", &comm_);
|
|
|
@@ -1503,21 +1654,35 @@ template <class Real> class Quadrature {
|
|
|
SCTL_ASSERT(U_near_sing.Dim() == U_direct.Dim());
|
|
|
Profile::Toc();
|
|
|
|
|
|
- if (U.Dim() != elements.NElem() * kernel.TrgDim()) {
|
|
|
- U.ReInit(elements.NElem() * kernel.TrgDim());
|
|
|
+ const Long dof = U_direct.Dim() / (elements.NElem() * PotentialBasis::Size() * kernel.TrgDim());
|
|
|
+ SCTL_ASSERT(U_direct .Dim() == elements.NElem() * PotentialBasis::Size() * dof * kernel.TrgDim());
|
|
|
+ SCTL_ASSERT(U_near_sing.Dim() == elements.NElem() * PotentialBasis::Size() * dof * kernel.TrgDim());
|
|
|
+
|
|
|
+ if (U.Dim() != elements.NElem() * dof * kernel.TrgDim()) {
|
|
|
+ U.ReInit(elements.NElem() * dof * kernel.TrgDim());
|
|
|
}
|
|
|
for (int i = 0; i < elements.NElem(); i++) {
|
|
|
for (int j = 0; j < PotentialBasis::Size(); j++) {
|
|
|
- for (int k = 0; k < kernel.TrgDim(); k++) {
|
|
|
- Real& U_ = U[i*kernel.TrgDim()+k][j];
|
|
|
+ for (int k = 0; k < dof*kernel.TrgDim(); k++) {
|
|
|
+ Real& U_ = U[i*dof*kernel.TrgDim()+k][j];
|
|
|
U_ = 0;
|
|
|
- U_ += U_direct [(i*PotentialBasis::Size()+j)*kernel.TrgDim()+k];
|
|
|
- U_ += U_near_sing[(i*PotentialBasis::Size()+j)*kernel.TrgDim()+k];
|
|
|
- U_ += U_singular[i*kernel.TrgDim()+k][j];
|
|
|
+ U_ += U_direct [(i*PotentialBasis::Size()+j)*dof*kernel.TrgDim()+k];
|
|
|
+ U_ += U_near_sing[(i*PotentialBasis::Size()+j)*dof*kernel.TrgDim()+k];
|
|
|
U_ *= kernel.template ScaleFactor<Real>();
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
+ if (U_singular.Dim(1)) {
|
|
|
+ SCTL_ASSERT(U_singular.Dim(0) == elements.NElem() * dof * kernel.TrgDim());
|
|
|
+ SCTL_ASSERT(U_singular.Dim(1) == PotentialBasis::Size());
|
|
|
+ for (int i = 0; i < elements.NElem(); i++) {
|
|
|
+ for (int j = 0; j < PotentialBasis::Size(); j++) {
|
|
|
+ for (int k = 0; k < dof*kernel.TrgDim(); k++) {
|
|
|
+ U[i*dof*kernel.TrgDim()+k][j] += U_singular[i*dof*kernel.TrgDim()+k][j] * kernel.template ScaleFactor<Real>();
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
Profile::Toc();
|
|
|
}
|
|
|
|
|
|
@@ -1539,6 +1704,11 @@ template <class Real> class Quadrature {
|
|
|
SCTL_ASSERT(U_near_sing.Dim() == U_direct.Dim());
|
|
|
Profile::Toc();
|
|
|
|
|
|
+
|
|
|
+ Long Nt = Xt_.Dim() / ElemList::CoordDim();
|
|
|
+ const Long dof = U_direct.Dim() / (Nt * kernel.TrgDim());
|
|
|
+ SCTL_ASSERT(U_direct.Dim() == Nt * dof * kernel.TrgDim());
|
|
|
+
|
|
|
if (U.Dim() != U_direct.Dim()) {
|
|
|
U.ReInit(U_direct.Dim());
|
|
|
}
|
|
|
@@ -1546,11 +1716,13 @@ template <class Real> class Quadrature {
|
|
|
U[i] = (U_direct[i] + U_near_sing[i]) * kernel.template ScaleFactor<Real>();
|
|
|
}
|
|
|
if (U_singular.Dim(1)) {
|
|
|
+ SCTL_ASSERT(U_singular.Dim(0) == elements.NElem() * dof * kernel.TrgDim());
|
|
|
+ const Long Nnodes = U_singular.Dim(1);
|
|
|
for (int i = 0; i < elements.NElem(); i++) {
|
|
|
- for (int j = 0; j < U_singular.Dim(1); j++) {
|
|
|
- for (int k = 0; k < kernel.TrgDim(); k++) {
|
|
|
- Real& U_ = U[(i*U_singular.Dim(1)+j)*kernel.TrgDim()+k];
|
|
|
- U_ += U_singular[i*kernel.TrgDim()+k][j] * kernel.template ScaleFactor<Real>();
|
|
|
+ for (int j = 0; j < Nnodes; j++) {
|
|
|
+ for (int k = 0; k < dof*kernel.TrgDim(); k++) {
|
|
|
+ Real& U_ = U[(i*Nnodes+j)*dof*kernel.TrgDim()+k];
|
|
|
+ U_ += U_singular[i*dof*kernel.TrgDim()+k][j] * kernel.template ScaleFactor<Real>();
|
|
|
}
|
|
|
}
|
|
|
}
|
|
|
@@ -1583,7 +1755,7 @@ template <class Real> class Quadrature {
|
|
|
for (long ii = start; ii < end; ii++) {
|
|
|
long i = ii / Np;
|
|
|
long j = ii % Np;
|
|
|
- for (int k = 0; k < nodes.Dim(1); k++) {
|
|
|
+ for (int k = 0; k < ElemList::CoordBasis::Size(); k++) {
|
|
|
Real X, Y, Z;
|
|
|
Real theta = 2 * const_pi<Real>() * (i + nodes[0][k]) / Nt;
|
|
|
Real phi = 2 * const_pi<Real>() * (j + nodes[1][k]) / Np;
|
|
|
@@ -1706,6 +1878,84 @@ template <class Real> class Quadrature {
|
|
|
Profile::print(&comm);
|
|
|
}
|
|
|
|
|
|
+ static void test1() {
|
|
|
+ const Comm& comm = Comm::World();
|
|
|
+ constexpr Integer ORDER = 15;
|
|
|
+ Integer order_singular = 20;
|
|
|
+ Integer order_direct = 20;
|
|
|
+
|
|
|
+ constexpr Integer COORD_DIM = 3;
|
|
|
+ constexpr Integer ELEM_DIM = COORD_DIM-1;
|
|
|
+ using ElemList = ElemList<COORD_DIM, Basis<Real, ELEM_DIM, ORDER>>;
|
|
|
+ using DensityBasis = Basis<Real, ELEM_DIM, ORDER>;
|
|
|
+ using PotentialBasis = Basis<Real, ELEM_DIM, ORDER>;
|
|
|
+
|
|
|
+ int np = comm.Size();
|
|
|
+ int rank = comm.Rank();
|
|
|
+ auto build_sphere = [rank,np](ElemList& elements, Real X, Real Y, Real Z, Real R){
|
|
|
+ auto nodes = ElemList::CoordBasis::Nodes();
|
|
|
+
|
|
|
+ long start = 2*COORD_DIM*(rank+0)/np;
|
|
|
+ long end = 2*COORD_DIM*(rank+1)/np;
|
|
|
+ elements.ReInit(end - start);
|
|
|
+ for (long ii = start; ii < end; ii++) {
|
|
|
+ long i = ii / 2;
|
|
|
+ long j = ii % 2;
|
|
|
+ for (int k = 0; k < ElemList::CoordBasis::Size(); k++) {
|
|
|
+ Real coord[COORD_DIM];
|
|
|
+ coord[(i+0)%COORD_DIM] = (j ? -1.0 : 1.0);
|
|
|
+ coord[(i+1)%COORD_DIM] = 2.0 * nodes[j?1:0][k] - 1.0;
|
|
|
+ coord[(i+2)%COORD_DIM] = 2.0 * nodes[j?0:1][k] - 1.0;
|
|
|
+ Real R0 = sqrt<Real>(coord[0]*coord[0] + coord[1]*coord[1] + coord[2]*coord[2]);
|
|
|
+
|
|
|
+ elements(ii-start,0)[k] = X + R * coord[0] / R0;
|
|
|
+ elements(ii-start,1)[k] = Y + R * coord[1] / R0;
|
|
|
+ elements(ii-start,2)[k] = Z + R * coord[2] / R0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ ElemList elements;
|
|
|
+ build_sphere(elements, 0.0, 0.0, 0.0, 1.00);
|
|
|
+
|
|
|
+ Vector<DensityBasis> density_sl;
|
|
|
+ { // Set density_sl
|
|
|
+ std::function<void(Real*,Real*,Real*)> sigma = [](Real* U, Real* X, Real* Xn) {
|
|
|
+ Real R = sqrt(X[0]*X[0]+X[1]*X[1]+X[2]*X[2]);
|
|
|
+ Real sinp = sqrt(X[1]*X[1] + X[2]*X[2]) / R;
|
|
|
+ Real cosp = -X[0] / R;
|
|
|
+
|
|
|
+ U[0] = -1.5;
|
|
|
+ U[1] = 0;
|
|
|
+ U[2] = 0;
|
|
|
+ };
|
|
|
+ DiscretizeSurfaceFn<COORD_DIM,3>(density_sl, elements, sigma);
|
|
|
+ }
|
|
|
+
|
|
|
+ GenericKernel<Stokes3D_DxU> Stokes_DxU;
|
|
|
+ GenericKernel<Stokes3D_FxU> Stokes_FxU;
|
|
|
+
|
|
|
+ Profile::Enable(true);
|
|
|
+ if (1) {
|
|
|
+ Vector<PotentialBasis> U;
|
|
|
+ Quadrature<Real> quadrature_FxU;
|
|
|
+ quadrature_FxU.Setup<DensityBasis, PotentialBasis>(elements, Stokes_FxU, order_singular, order_direct, -1.0, comm);
|
|
|
+ quadrature_FxU.Eval(U, elements, density_sl, Stokes_FxU);
|
|
|
+
|
|
|
+ { // Write VTK output
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(elements, U, ORDER);
|
|
|
+ vtu.WriteVTK("U", comm);
|
|
|
+ }
|
|
|
+ { // Write VTK output
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(elements, density_sl, ORDER);
|
|
|
+ vtu.WriteVTK("sigma", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Profile::print(&comm);
|
|
|
+ }
|
|
|
+
|
|
|
private:
|
|
|
|
|
|
static void scan(Vector<Long>& dsp, const Vector<Long>& cnt) {
|
|
|
@@ -1784,6 +2034,2296 @@ template <class Real> class Quadrature {
|
|
|
Comm comm_;
|
|
|
};
|
|
|
|
|
|
+
|
|
|
+template <class Real, Integer ORDER=10> class Stellarator {
|
|
|
+ private:
|
|
|
+ static constexpr Integer COORD_DIM = 3;
|
|
|
+ static constexpr Integer ELEM_DIM = COORD_DIM-1;
|
|
|
+ using ElemBasis = Basis<Real, ELEM_DIM, ORDER>;
|
|
|
+ using ElemLst = ElemList<COORD_DIM, ElemBasis>;
|
|
|
+
|
|
|
+ struct Laplace3D_dUxF {
|
|
|
+ template <class ValueType> static constexpr ValueType ScaleFactor() {
|
|
|
+ return 1 / (4 * const_pi<ValueType>());
|
|
|
+ }
|
|
|
+ template <class ValueType> static void Eval(ValueType (&u)[3][1], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) {
|
|
|
+ ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
|
|
|
+ ValueType rinv = (r2>1e-16 ? 1/sqrt<ValueType>(r2) : 0);
|
|
|
+ ValueType rinv3 = rinv * rinv * rinv;
|
|
|
+ u[0][0] = -r[0] * rinv3;
|
|
|
+ u[1][0] = -r[1] * rinv3;
|
|
|
+ u[2][0] = -r[2] * rinv3;
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ struct BiotSavart3D {
|
|
|
+ template <class ValueType> static constexpr ValueType ScaleFactor() {
|
|
|
+ return 1 / (4 * const_pi<ValueType>());
|
|
|
+ }
|
|
|
+ template <class ValueType> static void Eval(ValueType (&u)[3][3], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) {
|
|
|
+ ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
|
|
|
+ ValueType rinv = (r2>1e-16 ? 1/sqrt<ValueType>(r2) : 0);
|
|
|
+ ValueType rinv3 = rinv * rinv * rinv;
|
|
|
+ u[0][0] = (0) * rinv3; u[0][1] = -r[2] * rinv3; u[0][2] = r[1] * rinv3;
|
|
|
+ u[1][0] = r[2] * rinv3; u[1][1] = (0) * rinv3; u[1][2] = -r[0] * rinv3;
|
|
|
+ u[2][0] = -r[1] * rinv3; u[2][1] = r[0] * rinv3; u[2][2] = (0) * rinv3;
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ struct Laplace3D_dUxD {
|
|
|
+ template <class ValueType> static constexpr ValueType ScaleFactor() {
|
|
|
+ return 1 / (4 * const_pi<ValueType>());
|
|
|
+ }
|
|
|
+ template <class ValueType> static void Eval(ValueType (&u)[3][3], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) {
|
|
|
+ ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
|
|
|
+ ValueType rinv = (r2>1e-16 ? 1/sqrt<ValueType>(r2) : 0);
|
|
|
+ ValueType rdotn = r[0]*n[0] + r[1]*n[1] + r[2]*n[2];
|
|
|
+ ValueType rinv2 = rinv * rinv;
|
|
|
+ ValueType rinv3 = rinv * rinv2;
|
|
|
+ ValueType rinv5 = rinv3 * rinv2;
|
|
|
+
|
|
|
+ u[0][0] = -1 * rinv3 + 3 * r[0] * r[0] * rinv5;
|
|
|
+ u[0][1] = -0 * rinv3 + 3 * r[0] * r[1] * rinv5;
|
|
|
+ u[0][2] = -0 * rinv3 + 3 * r[0] * r[2] * rinv5;
|
|
|
+
|
|
|
+ u[1][0] = -0 * rinv3 + 3 * r[1] * r[0] * rinv5;
|
|
|
+ u[1][1] = -1 * rinv3 + 3 * r[1] * r[1] * rinv5;
|
|
|
+ u[1][2] = -0 * rinv3 + 3 * r[1] * r[2] * rinv5;
|
|
|
+
|
|
|
+ u[2][0] = -0 * rinv3 + 3 * r[2] * r[0] * rinv5;
|
|
|
+ u[2][1] = -0 * rinv3 + 3 * r[2] * r[1] * rinv5;
|
|
|
+ u[2][2] = -1 * rinv3 + 3 * r[2] * r[2] * rinv5;
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ struct Laplace3D_DxdU {
|
|
|
+ template <class ValueType> static constexpr ValueType ScaleFactor() {
|
|
|
+ return 1 / (4 * const_pi<ValueType>());
|
|
|
+ }
|
|
|
+ template <class ValueType> static void Eval(ValueType (&u)[1][3], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) {
|
|
|
+ ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
|
|
|
+ ValueType rinv = (r2>1e-16 ? 1/sqrt<ValueType>(r2) : 0);
|
|
|
+ ValueType rdotn = r[0]*n[0] + r[1]*n[1] + r[2]*n[2];
|
|
|
+ ValueType rinv2 = rinv * rinv;
|
|
|
+ ValueType rinv3 = rinv * rinv2;
|
|
|
+ ValueType rinv5 = rinv3 * rinv2;
|
|
|
+ u[0][0] = -n[0] * rinv3 + 3*rdotn * r[0] * rinv5;
|
|
|
+ u[0][1] = -n[1] * rinv3 + 3*rdotn * r[1] * rinv5;
|
|
|
+ u[0][2] = -n[2] * rinv3 + 3*rdotn * r[2] * rinv5;
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ struct Laplace3D_Fxd2U {
|
|
|
+ template <class ValueType> static constexpr ValueType ScaleFactor() {
|
|
|
+ return 1 / (4 * const_pi<ValueType>());
|
|
|
+ }
|
|
|
+ template <class ValueType> static void Eval(ValueType (&u)[1][9], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) {
|
|
|
+ ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
|
|
|
+ ValueType rinv = (r2>1e-16 ? 1/sqrt<ValueType>(r2) : 0);
|
|
|
+ ValueType rinv2 = rinv * rinv;
|
|
|
+ ValueType rinv3 = rinv * rinv2;
|
|
|
+ ValueType rinv5 = rinv3 * rinv2;
|
|
|
+
|
|
|
+ u[0][0+3*0] = -1 * rinv3 + 3 * r[0] * r[0] * rinv5;
|
|
|
+ u[0][1+3*0] = -0 * rinv3 + 3 * r[0] * r[1] * rinv5;
|
|
|
+ u[0][2+3*0] = -0 * rinv3 + 3 * r[0] * r[2] * rinv5;
|
|
|
+
|
|
|
+ u[0][0+3*1] = -0 * rinv3 + 3 * r[1] * r[0] * rinv5;
|
|
|
+ u[0][1+3*1] = -1 * rinv3 + 3 * r[1] * r[1] * rinv5;
|
|
|
+ u[0][2+3*1] = -0 * rinv3 + 3 * r[1] * r[2] * rinv5;
|
|
|
+
|
|
|
+ u[0][0+3*2] = -0 * rinv3 + 3 * r[2] * r[0] * rinv5;
|
|
|
+ u[0][1+3*2] = -0 * rinv3 + 3 * r[2] * r[1] * rinv5;
|
|
|
+ u[0][2+3*2] = -1 * rinv3 + 3 * r[2] * r[2] * rinv5;
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ public:
|
|
|
+ Stellarator(const Vector<Long>& NtNp = Vector<Long>()) {
|
|
|
+ NtNp_ = NtNp;
|
|
|
+ Long Nsurf = NtNp_.Dim() / 2;
|
|
|
+ SCTL_ASSERT(Nsurf*2 == NtNp_.Dim());
|
|
|
+
|
|
|
+ Long Nelem = 0;
|
|
|
+ elem_dsp.ReInit(Nsurf);
|
|
|
+ if (elem_dsp.Dim()) elem_dsp[0] = 0;
|
|
|
+ for (Long i = 0; i < Nsurf; i++) {
|
|
|
+ Nelem += NtNp_[i*2+0]*NtNp_[i*2+1];
|
|
|
+ if (i+1 < Nsurf) elem_dsp[i+1] = elem_dsp[i] + NtNp_[i*2+0]*NtNp_[i*2+1];
|
|
|
+ }
|
|
|
+ elements.ReInit(Nelem);
|
|
|
+ for (Long i = 0; i < Nsurf; i++) {
|
|
|
+ InitSurf(i);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Long ElemIdx(Long s, Long t, Long p) {
|
|
|
+ SCTL_ASSERT(0 <= s && s < elem_dsp.Dim());
|
|
|
+ SCTL_ASSERT(0 <= t && t < NtNp_[s*2+0]);
|
|
|
+ SCTL_ASSERT(0 <= p && p < NtNp_[s*2+1]);
|
|
|
+ return elem_dsp[s] + t*NtNp_[s*2+1] + p;
|
|
|
+ }
|
|
|
+ ElemBasis& Elem(Long elem, Integer dim) {
|
|
|
+ return elements(elem,dim);
|
|
|
+ }
|
|
|
+ const ElemBasis& Elem(Long elem, Integer dim) const {
|
|
|
+ return elements(elem,dim);
|
|
|
+ }
|
|
|
+ const ElemLst& GetElemList() {
|
|
|
+ return elements;
|
|
|
+ }
|
|
|
+
|
|
|
+ static void test_() {
|
|
|
+ constexpr Integer order_singular = 20;
|
|
|
+ constexpr Integer order_direct = 35;
|
|
|
+ Comm comm = Comm::World();
|
|
|
+ Profile::Enable(true);
|
|
|
+
|
|
|
+ Stellarator<Real,ORDER> S;
|
|
|
+ { // Set S
|
|
|
+ Vector<Real> X(COORD_DIM);
|
|
|
+ Vector<Real> R(1);
|
|
|
+ X = 0;
|
|
|
+ R = 1;
|
|
|
+ SCTL_ASSERT(X.Dim() == R.Dim() * COORD_DIM);
|
|
|
+ Long N = R.Dim();
|
|
|
+ S.elements.ReInit(2*COORD_DIM*N);
|
|
|
+ auto nodes = ElemLst::CoordBasis::Nodes();
|
|
|
+ for (Long l = 0; l < N; l++) {
|
|
|
+ for (Integer i = 0; i < COORD_DIM; i++) {
|
|
|
+ for (Integer j = 0; j < 2; j++) {
|
|
|
+ for (int k = 0; k < ElemLst::CoordBasis::Size(); k++) {
|
|
|
+ Real coord[COORD_DIM];
|
|
|
+ coord[(i+0)%COORD_DIM] = (j ? -1.0 : 1.0);
|
|
|
+ coord[(i+1)%COORD_DIM] = 2.0 * nodes[j?1:0][k] - 1.0;
|
|
|
+ coord[(i+2)%COORD_DIM] = 2.0 * nodes[j?0:1][k] - 1.0;
|
|
|
+ Real R0 = sqrt<Real>(coord[0]*coord[0] + coord[1]*coord[1] + coord[2]*coord[2]);
|
|
|
+
|
|
|
+ S.elements((l*COORD_DIM+i)*2+j,0)[k] = X[l*COORD_DIM+0] + R[l] * coord[0] / R0;
|
|
|
+ S.elements((l*COORD_DIM+i)*2+j,1)[k] = X[l*COORD_DIM+1] + R[l] * coord[1] / R0;
|
|
|
+ S.elements((l*COORD_DIM+i)*2+j,2)[k] = X[l*COORD_DIM+2] + R[l] * coord[2] / R0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ S.elem_dsp.ReInit(1);
|
|
|
+ S.elem_dsp = 0;
|
|
|
+ }
|
|
|
+
|
|
|
+ S.quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
|
|
|
+ //S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
|
|
|
+
|
|
|
+ { // test Fxd2U
|
|
|
+ Vector<ElemBasis> U, sigma(S.elements.NElem());
|
|
|
+ sigma = 1;
|
|
|
+ sigma[0] = 1;
|
|
|
+ S.quadrature_Fxd2U.Eval(U, S.GetElemList(), sigma, S.Laplace_Fxd2U);
|
|
|
+ //S.quadrature_FxdU.Eval(U, S.GetElemList(), sigma, S.Laplace_FxdU);
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), U, ORDER);
|
|
|
+ vtu.WriteVTK("test", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Profile::print(&comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ static void test() {
|
|
|
+ constexpr Integer order_singular = 15;
|
|
|
+ constexpr Integer order_direct = 35;
|
|
|
+ Comm comm = Comm::World();
|
|
|
+ Profile::Enable(true);
|
|
|
+
|
|
|
+ Stellarator<Real,ORDER> S;
|
|
|
+ { // Init S
|
|
|
+ Vector<Long> NtNp;
|
|
|
+ NtNp.PushBack(20);
|
|
|
+ NtNp.PushBack(4);
|
|
|
+ S = Stellarator<Real,ORDER>(NtNp);
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> normal, area_elem;
|
|
|
+ auto compute_dot_prod = [](const Vector<ElemBasis>& A, const Vector<ElemBasis>& B) {
|
|
|
+ const Long Nelem = A.Dim() / COORD_DIM;
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ SCTL_ASSERT(A.Dim() == Nelem * COORD_DIM);
|
|
|
+ SCTL_ASSERT(B.Dim() == Nelem * COORD_DIM);
|
|
|
+ Vector<ElemBasis> AdotB(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Real a_dot_b = 0;
|
|
|
+ a_dot_b += A[i*COORD_DIM+0][j]*B[i*COORD_DIM+0][j];
|
|
|
+ a_dot_b += A[i*COORD_DIM+1][j]*B[i*COORD_DIM+1][j];
|
|
|
+ a_dot_b += A[i*COORD_DIM+2][j]*B[i*COORD_DIM+2][j];
|
|
|
+ AdotB[i][j] = a_dot_b;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return AdotB;
|
|
|
+ };
|
|
|
+ auto compute_inner_prod = [&S, &area_elem](const Vector<ElemBasis>& A, const Vector<ElemBasis>& B) {
|
|
|
+ const auto& quad_wts = ElemBasis::QuadWts();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ const Long dof = B.Dim() / Nelem;
|
|
|
+ Real sum = 0;
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Real AdotB = 0;
|
|
|
+ for (Long k = 0; k < dof; k++) {
|
|
|
+ AdotB += A[i*dof+k][j] * B[i*dof+k][j];
|
|
|
+ }
|
|
|
+ sum += AdotB * area_elem[i][j] * quad_wts[j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return sum;
|
|
|
+ };
|
|
|
+ auto compute_norm_area_elem = [&S](Vector<ElemBasis>& normal, Vector<ElemBasis>& area_elem){ // Set normal, area_elem
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+
|
|
|
+ area_elem.ReInit(Nelem);
|
|
|
+ normal.ReInit(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM> x, n;
|
|
|
+ Tensor<Real,true,COORD_DIM,2> dx;
|
|
|
+ x(0) = X[i*COORD_DIM+0][j];
|
|
|
+ x(1) = X[i*COORD_DIM+1][j];
|
|
|
+ x(2) = X[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dx(0,0) = dX[i*COORD_DIM*2+0][j];
|
|
|
+ dx(0,1) = dX[i*COORD_DIM*2+1][j];
|
|
|
+ dx(1,0) = dX[i*COORD_DIM*2+2][j];
|
|
|
+ dx(1,1) = dX[i*COORD_DIM*2+3][j];
|
|
|
+ dx(2,0) = dX[i*COORD_DIM*2+4][j];
|
|
|
+ dx(2,1) = dX[i*COORD_DIM*2+5][j];
|
|
|
+
|
|
|
+ n(0) = dx(1,0) * dx(2,1) - dx(2,0) * dx(1,1);
|
|
|
+ n(1) = dx(2,0) * dx(0,1) - dx(0,0) * dx(2,1);
|
|
|
+ n(2) = dx(0,0) * dx(1,1) - dx(1,0) * dx(0,1);
|
|
|
+ Real area_elem_ = sqrt<Real>(n(0)*n(0) + n(1)*n(1) + n(2)*n(2));
|
|
|
+ Real ooae = 1 / area_elem_;
|
|
|
+ n(0) *= ooae;
|
|
|
+ n(1) *= ooae;
|
|
|
+ n(2) *= ooae;
|
|
|
+
|
|
|
+ normal[i*COORD_DIM+0][j] = n(0);
|
|
|
+ normal[i*COORD_DIM+1][j] = n(1);
|
|
|
+ normal[i*COORD_DIM+2][j] = n(2);
|
|
|
+ area_elem[i][j] = area_elem_;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ auto compute_poloidal_circulation = [&S] (const Vector<ElemBasis>& B) {
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ const auto& quad_wts = Basis<Real,1,ORDER>::QuadWts();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+ const Long Nt = 40;
|
|
|
+ const Long Np = 8;
|
|
|
+ for (Long t = 0; t < Nt; t++) {
|
|
|
+ for (Long j = 0; j < ORDER; j++) {
|
|
|
+ Real sum = 0;
|
|
|
+ for (Long p = 0; p < Np; p++) {
|
|
|
+ for (Long i = 0; i < ORDER; i++) {
|
|
|
+ Long elem_idx = t*Np+p;
|
|
|
+ Long node_idx = i*ORDER+j;
|
|
|
+
|
|
|
+ Tensor<Real,true,COORD_DIM,2> dx;
|
|
|
+ dx(0,0) = dX[elem_idx*COORD_DIM*2+0][node_idx];
|
|
|
+ dx(0,1) = dX[elem_idx*COORD_DIM*2+1][node_idx];
|
|
|
+ dx(1,0) = dX[elem_idx*COORD_DIM*2+2][node_idx];
|
|
|
+ dx(1,1) = dX[elem_idx*COORD_DIM*2+3][node_idx];
|
|
|
+ dx(2,0) = dX[elem_idx*COORD_DIM*2+4][node_idx];
|
|
|
+ dx(2,1) = dX[elem_idx*COORD_DIM*2+5][node_idx];
|
|
|
+
|
|
|
+ Tensor<Real,true,COORD_DIM> b;
|
|
|
+ b(0) = B[elem_idx*COORD_DIM+0][node_idx];
|
|
|
+ b(1) = B[elem_idx*COORD_DIM+1][node_idx];
|
|
|
+ b(2) = B[elem_idx*COORD_DIM+2][node_idx];
|
|
|
+
|
|
|
+ sum += (b(0)*dx(0,1) + b(1)*dx(1,1) + b(2)*dx(2,1)) * quad_wts[i];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ std::cout<<sum<<' ';
|
|
|
+ }
|
|
|
+ }
|
|
|
+ std::cout<<'\n';
|
|
|
+ };
|
|
|
+ auto compute_toroidal_circulation = [&S] (const Vector<ElemBasis>& B) {
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ const auto& quad_wts = Basis<Real,1,ORDER>::QuadWts();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+ const Long Nt = 40;
|
|
|
+ const Long Np = 8;
|
|
|
+ for (Long p = 0; p < Np; p++) {
|
|
|
+ for (Long i = 0; i < ORDER; i++) {
|
|
|
+ Real sum = 0;
|
|
|
+ for (Long t = 0; t < Nt; t++) {
|
|
|
+ for (Long j = 0; j < ORDER; j++) {
|
|
|
+ Long elem_idx = t*Np+p;
|
|
|
+ Long node_idx = i*ORDER+j;
|
|
|
+
|
|
|
+ Tensor<Real,true,COORD_DIM,2> dx;
|
|
|
+ dx(0,0) = dX[elem_idx*COORD_DIM*2+0][node_idx];
|
|
|
+ dx(0,1) = dX[elem_idx*COORD_DIM*2+1][node_idx];
|
|
|
+ dx(1,0) = dX[elem_idx*COORD_DIM*2+2][node_idx];
|
|
|
+ dx(1,1) = dX[elem_idx*COORD_DIM*2+3][node_idx];
|
|
|
+ dx(2,0) = dX[elem_idx*COORD_DIM*2+4][node_idx];
|
|
|
+ dx(2,1) = dX[elem_idx*COORD_DIM*2+5][node_idx];
|
|
|
+
|
|
|
+ Tensor<Real,true,COORD_DIM> b;
|
|
|
+ b(0) = B[elem_idx*COORD_DIM+0][node_idx];
|
|
|
+ b(1) = B[elem_idx*COORD_DIM+1][node_idx];
|
|
|
+ b(2) = B[elem_idx*COORD_DIM+2][node_idx];
|
|
|
+
|
|
|
+ sum += (b(0)*dx(0,0) + b(1)*dx(1,0) + b(2)*dx(2,0)) * quad_wts[j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ std::cout<<sum<<' ';
|
|
|
+ }
|
|
|
+ }
|
|
|
+ std::cout<<'\n';
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_poloidal_circulation_ = [&S,&area_elem] (const Vector<ElemBasis>& B) {
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ const auto& quad_wts = ElemBasis::QuadWts();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+ Real sum = 0;
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM,2> dx;
|
|
|
+ dx(0,0) = dX[i*COORD_DIM*2+0][j];
|
|
|
+ dx(0,1) = dX[i*COORD_DIM*2+1][j];
|
|
|
+ dx(1,0) = dX[i*COORD_DIM*2+2][j];
|
|
|
+ dx(1,1) = dX[i*COORD_DIM*2+3][j];
|
|
|
+ dx(2,0) = dX[i*COORD_DIM*2+4][j];
|
|
|
+ dx(2,1) = dX[i*COORD_DIM*2+5][j];
|
|
|
+
|
|
|
+ Tensor<Real,true,COORD_DIM> b;
|
|
|
+ b(0) = B[i*COORD_DIM+0][j];
|
|
|
+ b(1) = B[i*COORD_DIM+1][j];
|
|
|
+ b(2) = B[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ Real s = 1/area_elem[i][j];
|
|
|
+ sum += (b(0)*dx(0,1) + b(1)*dx(1,1) + b(2)*dx(2,1)) * s * area_elem[i][j] * quad_wts[j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return sum;
|
|
|
+ };
|
|
|
+ auto compute_toroidal_circulation_ = [&S,&area_elem] (const Vector<ElemBasis>& B) {
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ const auto& quad_wts = ElemBasis::QuadWts();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+ Real sum = 0;
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM,2> dx;
|
|
|
+ dx(0,0) = dX[i*COORD_DIM*2+0][j];
|
|
|
+ dx(0,1) = dX[i*COORD_DIM*2+1][j];
|
|
|
+ dx(1,0) = dX[i*COORD_DIM*2+2][j];
|
|
|
+ dx(1,1) = dX[i*COORD_DIM*2+3][j];
|
|
|
+ dx(2,0) = dX[i*COORD_DIM*2+4][j];
|
|
|
+ dx(2,1) = dX[i*COORD_DIM*2+5][j];
|
|
|
+
|
|
|
+ Tensor<Real,true,COORD_DIM> b;
|
|
|
+ b(0) = B[i*COORD_DIM+0][j];
|
|
|
+ b(1) = B[i*COORD_DIM+1][j];
|
|
|
+ b(2) = B[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ Real s = 1/area_elem[i][j];
|
|
|
+ sum += (b(0)*dx(0,0) + b(1)*dx(1,0) + b(2)*dx(2,0)) * s * area_elem[i][j] * quad_wts[j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return sum;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_grad_adj = [&S,&area_elem] (const Vector<ElemBasis>& V) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> du_dX(Nelem*COORD_DIM*2);
|
|
|
+ { // Set du_dX
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, S.GetElemList().ElemVector());
|
|
|
+
|
|
|
+ auto inv2x2 = [](Tensor<Real, true, 2, 2> M) {
|
|
|
+ Tensor<Real, true, 2, 2> Mout;
|
|
|
+ Real oodet = 1 / (M(0,0) * M(1,1) - M(0,1) * M(1,0));
|
|
|
+ Mout(0,0) = M(1,1) * oodet;
|
|
|
+ Mout(0,1) = -M(0,1) * oodet;
|
|
|
+ Mout(1,0) = -M(1,0) * oodet;
|
|
|
+ Mout(1,1) = M(0,0) * oodet;
|
|
|
+ return Mout;
|
|
|
+ };
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real, true, 3, 2> dX_du;
|
|
|
+ dX_du(0,0) = dX[(i*COORD_DIM+0)*2+0][j];
|
|
|
+ dX_du(1,0) = dX[(i*COORD_DIM+1)*2+0][j];
|
|
|
+ dX_du(2,0) = dX[(i*COORD_DIM+2)*2+0][j];
|
|
|
+ dX_du(0,1) = dX[(i*COORD_DIM+0)*2+1][j];
|
|
|
+ dX_du(1,1) = dX[(i*COORD_DIM+1)*2+1][j];
|
|
|
+ dX_du(2,1) = dX[(i*COORD_DIM+2)*2+1][j];
|
|
|
+
|
|
|
+ Tensor<Real, true, 2, 2> G; // = dX_du.Transpose() * dX_du;
|
|
|
+ G(0,0) = dX_du(0,0) * dX_du(0,0) + dX_du(1,0) * dX_du(1,0) + dX_du(2,0) * dX_du(2,0);
|
|
|
+ G(0,1) = dX_du(0,0) * dX_du(0,1) + dX_du(1,0) * dX_du(1,1) + dX_du(2,0) * dX_du(2,1);
|
|
|
+ G(1,0) = dX_du(0,1) * dX_du(0,0) + dX_du(1,1) * dX_du(1,0) + dX_du(2,1) * dX_du(2,0);
|
|
|
+ G(1,1) = dX_du(0,1) * dX_du(0,1) + dX_du(1,1) * dX_du(1,1) + dX_du(2,1) * dX_du(2,1);
|
|
|
+
|
|
|
+ Tensor<Real, true, 2, 2> Ginv = inv2x2(G);
|
|
|
+ du_dX[(i*COORD_DIM+0)*2+0][j] = Ginv(0,0) * dX_du(0,0) + Ginv(0,1) * dX_du(0,1);
|
|
|
+ du_dX[(i*COORD_DIM+1)*2+0][j] = Ginv(0,0) * dX_du(1,0) + Ginv(0,1) * dX_du(1,1);
|
|
|
+ du_dX[(i*COORD_DIM+2)*2+0][j] = Ginv(0,0) * dX_du(2,0) + Ginv(0,1) * dX_du(2,1);
|
|
|
+ du_dX[(i*COORD_DIM+0)*2+1][j] = Ginv(1,0) * dX_du(0,0) + Ginv(1,1) * dX_du(0,1);
|
|
|
+ du_dX[(i*COORD_DIM+1)*2+1][j] = Ginv(1,0) * dX_du(1,0) + Ginv(1,1) * dX_du(1,1);
|
|
|
+ du_dX[(i*COORD_DIM+2)*2+1][j] = Ginv(1,0) * dX_du(2,0) + Ginv(1,1) * dX_du(2,1);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> dudX_V(Nelem*2);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dudX_V[i*2+0][j] = 0;
|
|
|
+ dudX_V[i*2+1][j] = 0;
|
|
|
+
|
|
|
+ dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+0)*2+0][j] * V[i*COORD_DIM+0][j];
|
|
|
+ dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+1)*2+0][j] * V[i*COORD_DIM+1][j];
|
|
|
+ dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+2)*2+0][j] * V[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+0)*2+1][j] * V[i*COORD_DIM+0][j];
|
|
|
+ dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+1)*2+1][j] * V[i*COORD_DIM+1][j];
|
|
|
+ dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+2)*2+1][j] * V[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> eye(Nnodes), Mgrad;
|
|
|
+ eye = 0;
|
|
|
+ for (Long i = 0; i < Nnodes; i++) eye[i][i] = 1;
|
|
|
+ ElemBasis::Grad(Mgrad, eye);
|
|
|
+
|
|
|
+ Vector<ElemBasis> grad_adj_V(Nelem);
|
|
|
+ const auto& quad_wts = ElemBasis::QuadWts();
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Real sum = 0;
|
|
|
+ for (Long k = 0; k < Nnodes; k++) {
|
|
|
+ sum += Mgrad[j*2+0][k] * dudX_V[i*2+0][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
|
|
|
+ sum += Mgrad[j*2+1][k] * dudX_V[i*2+1][k] * (area_elem[i][k] * quad_wts[k]) / (quad_wts[j] * area_elem[i][j]);
|
|
|
+ }
|
|
|
+ grad_adj_V[i][j] = -sum;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return grad_adj_V;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_B0 = [&S](const Real alpha) { // alpha/|r| \hat{\theta}
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> B0(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM> x, b0, axis;
|
|
|
+ x(0) = X[i*COORD_DIM+0][j];
|
|
|
+ x(1) = X[i*COORD_DIM+1][j];
|
|
|
+ x(2) = X[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ axis(0) = 0;
|
|
|
+ axis(1) = 0;
|
|
|
+ axis(2) = 1;
|
|
|
+ b0(0) = axis(1) * x(2) - axis(2) * x(1);
|
|
|
+ b0(1) = axis(2) * x(0) - axis(0) * x(2);
|
|
|
+ b0(2) = axis(0) * x(1) - axis(1) * x(0);
|
|
|
+ Real scale = 1 / (b0(0)*b0(0) + b0(1)*b0(1) + b0(2)*b0(2));
|
|
|
+ b0(0) *= scale;
|
|
|
+ b0(1) *= scale;
|
|
|
+ b0(2) *= scale;
|
|
|
+
|
|
|
+ B0[i*COORD_DIM+0][j] = alpha * b0(0);
|
|
|
+ B0[i*COORD_DIM+1][j] = alpha * b0(1);
|
|
|
+ B0[i*COORD_DIM+2][j] = alpha * b0(2);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return B0;
|
|
|
+ };
|
|
|
+ auto compute_dB0 = [&S](const Real alpha) { // alpha/|r| \hat{\theta}
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dB0(Nelem * COORD_DIM * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM> x;
|
|
|
+ x(0) = X[i*COORD_DIM+0][j];
|
|
|
+ x(1) = X[i*COORD_DIM+1][j];
|
|
|
+ x(2) = X[i*COORD_DIM+2][j];
|
|
|
+ Real R2inv = 1 / (x(0)*x(0) + x(1)*x(1));
|
|
|
+
|
|
|
+ dB0[(i*COORD_DIM+0)*COORD_DIM+0][j] = alpha * (2*x(0)*x(1) * R2inv*R2inv);
|
|
|
+ dB0[(i*COORD_DIM+0)*COORD_DIM+1][j] = alpha * (-R2inv + 2*x(1)*x(1) * R2inv*R2inv);
|
|
|
+ dB0[(i*COORD_DIM+0)*COORD_DIM+2][j] = 0;
|
|
|
+
|
|
|
+ dB0[(i*COORD_DIM+1)*COORD_DIM+0][j] = alpha * (R2inv - 2*x(0)*x(0) * R2inv*R2inv);
|
|
|
+ dB0[(i*COORD_DIM+1)*COORD_DIM+1][j] = alpha * (-2*x(0)*x(1) * R2inv*R2inv);
|
|
|
+ dB0[(i*COORD_DIM+1)*COORD_DIM+2][j] = 0;
|
|
|
+
|
|
|
+ dB0[(i*COORD_DIM+2)*COORD_DIM+0][j] = 0;
|
|
|
+ dB0[(i*COORD_DIM+2)*COORD_DIM+1][j] = 0;
|
|
|
+ dB0[(i*COORD_DIM+2)*COORD_DIM+2][j] = 0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return dB0;
|
|
|
+ };
|
|
|
+ auto compute_half_n_plus_dG = [&S, &normal](const Vector<ElemBasis>& sigma) { // B = n sigma/2 + dG[sigma]
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> B;
|
|
|
+ S.quadrature_FxdU.Eval(B, S.GetElemList(), sigma, S.Laplace_FxdU);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ for (Long k = 0; k < COORD_DIM; k++) {
|
|
|
+ B[i*COORD_DIM+k][j] -= 0.5*sigma[i][j]*normal[i*COORD_DIM+k][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return B;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_A21adj = [&S,&area_elem,&normal](Vector<Real>& A21adj_flux, Real flux) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> density(Nelem * COORD_DIM);
|
|
|
+ { // Set density
|
|
|
+ Vector<ElemBasis> dX;
|
|
|
+ ElemBasis::Grad(dX, S.GetElemList().ElemVector());
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM,2> dx;
|
|
|
+ dx(0,0) = dX[i*COORD_DIM*2+0][j];
|
|
|
+ dx(0,1) = dX[i*COORD_DIM*2+1][j];
|
|
|
+ dx(1,0) = dX[i*COORD_DIM*2+2][j];
|
|
|
+ dx(1,1) = dX[i*COORD_DIM*2+3][j];
|
|
|
+ dx(2,0) = dX[i*COORD_DIM*2+4][j];
|
|
|
+ dx(2,1) = dX[i*COORD_DIM*2+5][j];
|
|
|
+
|
|
|
+ Real s = 1 / (area_elem[i][j] * S.NtNp_[0]);
|
|
|
+ for (Long k = 0; k < COORD_DIM; k++) {
|
|
|
+ density[i*COORD_DIM+k][j] = dx(k,1) * s;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> Gdensity;
|
|
|
+ S.quadrature_FxU.Eval(Gdensity, S.GetElemList(), density, S.Laplace_FxU);
|
|
|
+
|
|
|
+ Vector<ElemBasis> nxGdensity(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) { // Set nxGdensity
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM> Gdensity_, n;
|
|
|
+ Gdensity_(0) = Gdensity[i*COORD_DIM+0][j];
|
|
|
+ Gdensity_(1) = Gdensity[i*COORD_DIM+1][j];
|
|
|
+ Gdensity_(2) = Gdensity[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ n(0) = normal[i*COORD_DIM+0][j];
|
|
|
+ n(1) = normal[i*COORD_DIM+1][j];
|
|
|
+ n(2) = normal[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ nxGdensity[i*COORD_DIM+0][j] = n(1) * Gdensity_(2) - n(2) * Gdensity_(1);
|
|
|
+ nxGdensity[i*COORD_DIM+1][j] = n(2) * Gdensity_(0) - n(0) * Gdensity_(2);
|
|
|
+ nxGdensity[i*COORD_DIM+2][j] = n(0) * Gdensity_(1) - n(1) * Gdensity_(0);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ S.quadrature_dUxF.Eval(A21adj_flux, S.GetElemList(), nxGdensity, S.Laplace_dUxF);
|
|
|
+ A21adj_flux *= flux;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_A11 = [&S,&normal,&compute_half_n_plus_dG,&compute_dot_prod](Vector<Real>& B_dot_n, const Vector<Real>& sigma) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ B_dot_n.ReInit(Nelem * Nnodes);
|
|
|
+ Vector<ElemBasis> sigma_(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i][j] = sigma[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Vector<ElemBasis> B_dot_n_ = compute_dot_prod(normal, compute_half_n_plus_dG(sigma_));
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ B_dot_n[i*Nnodes+j] = B_dot_n_[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ auto compute_A12 = [&S,&normal,&compute_dot_prod,&compute_B0](Vector<Real>& B_dot_n, const Real alpha) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ B_dot_n.ReInit(Nelem * Nnodes);
|
|
|
+ Vector<ElemBasis> B_dot_n_ = compute_dot_prod(normal, compute_B0(alpha));
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ B_dot_n[i*Nnodes+j] = B_dot_n_[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ auto compute_A21 = [&S,&normal,&compute_half_n_plus_dG,&compute_poloidal_circulation_,&compute_A21adj,&compute_inner_prod](const Vector<Real>& sigma) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> sigma_(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i][j] = sigma[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ if (0) {
|
|
|
+ Vector<ElemBasis> A21_(Nelem);
|
|
|
+ Vector<Real> A21(Nelem*Nnodes);
|
|
|
+ compute_A21adj(A21, 1);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ A21_[i][j] = A21[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return compute_inner_prod(A21_, sigma_);
|
|
|
+ } else {
|
|
|
+ Vector<ElemBasis> B = compute_half_n_plus_dG(sigma_);
|
|
|
+
|
|
|
+ Vector<ElemBasis> J(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) { // Set J
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM> b, n;
|
|
|
+ b(0) = B[i*COORD_DIM+0][j];
|
|
|
+ b(1) = B[i*COORD_DIM+1][j];
|
|
|
+ b(2) = B[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ n(0) = normal[i*COORD_DIM+0][j];
|
|
|
+ n(1) = normal[i*COORD_DIM+1][j];
|
|
|
+ n(2) = normal[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
|
|
|
+ J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
|
|
|
+ J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> A;
|
|
|
+ S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
|
|
|
+ return compute_poloidal_circulation_(A)/S.NtNp_[0];
|
|
|
+ }
|
|
|
+ };
|
|
|
+ auto compute_A22 = [&S,&compute_B0,&normal,&compute_poloidal_circulation_](const Real alpha) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> B = compute_B0(alpha);
|
|
|
+
|
|
|
+ Vector<ElemBasis> J(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) { // Set J
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,COORD_DIM> b, n;
|
|
|
+ b(0) = B[i*COORD_DIM+0][j];
|
|
|
+ b(1) = B[i*COORD_DIM+1][j];
|
|
|
+ b(2) = B[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ n(0) = normal[i*COORD_DIM+0][j];
|
|
|
+ n(1) = normal[i*COORD_DIM+1][j];
|
|
|
+ n(2) = normal[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1);
|
|
|
+ J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2);
|
|
|
+ J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> A;
|
|
|
+ S.quadrature_FxU.Eval(A, S.GetElemList(), J, S.Laplace_FxU);
|
|
|
+ return compute_poloidal_circulation_(A)/S.NtNp_[0];
|
|
|
+ };
|
|
|
+ auto compute_A = [&compute_A11,&compute_A12,&compute_A21,&compute_A22] (const Vector<Real>& x) {
|
|
|
+ const Vector<Real> sigma(x.Dim()-1,(Iterator<Real>)x.begin(),false);
|
|
|
+ const Real& alpha = x[x.Dim()-1];
|
|
|
+
|
|
|
+ Vector<Real> Ax;
|
|
|
+ Ax.ReInit(x.Dim());
|
|
|
+ Vector<Real> Bdotn(x.Dim()-1,Ax.begin(),false);
|
|
|
+ Real& flux = Ax[x.Dim()-1];
|
|
|
+
|
|
|
+ Vector<Real> Adotn_0, Adotn_1;
|
|
|
+ compute_A11(Adotn_0, sigma);
|
|
|
+ compute_A12(Adotn_1, alpha);
|
|
|
+ Bdotn = Adotn_0 + Adotn_1;
|
|
|
+
|
|
|
+ flux = compute_A21(sigma) + compute_A22(alpha);
|
|
|
+ return Ax;
|
|
|
+ };
|
|
|
+ auto compute_invA = [&S,&comm,&compute_A] (Vector<ElemBasis>& sigma, Real& alpha, Real flux) {
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&compute_A](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
|
|
|
+ (*Ax) = compute_A(x);
|
|
|
+ };
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<Real> rhs_(Nelem * Nnodes + 1);
|
|
|
+ rhs_ = 0;
|
|
|
+ rhs_[Nelem * Nnodes] = flux;
|
|
|
+
|
|
|
+ Vector<Real> x_(Nelem * Nnodes + 1);
|
|
|
+ x_ = 0;
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&x_, BIOp, rhs_, 1e-8, 50);
|
|
|
+
|
|
|
+ sigma.ReInit(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma[i][j] = x_[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ alpha = x_[Nelem * Nnodes];
|
|
|
+ };
|
|
|
+ auto compute_invA_ = [&S,&comm,&compute_A] (Vector<Real>& b) {
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&compute_A](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
|
|
|
+ (*Ax) = compute_A(x);
|
|
|
+ };
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<Real> x(b.Dim());
|
|
|
+ x = 0;
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&x, BIOp, b, 1e-8, 50);
|
|
|
+ return x;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_A11adj = [&S](Vector<Real>& U, const Vector<Real>& sigma) { // A11adj = I/2 + D
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> sigma_(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i][j] = sigma[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ S.quadrature_DxU.Eval(U, S.GetElemList(), sigma_, S.Laplace_DxU);
|
|
|
+ U = sigma*(-0.5) + U;
|
|
|
+ };
|
|
|
+ auto compute_A12adj = [&S,&compute_A12,&compute_inner_prod](const Vector<Real>& sigma_) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<Real> A12_sigma_;
|
|
|
+ compute_A12(A12_sigma_, 1);
|
|
|
+
|
|
|
+ Vector<ElemBasis> A12_sigma(Nelem), sigma(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma[i][j] = sigma_[i*Nnodes+j];
|
|
|
+ A12_sigma[i][j] = A12_sigma_[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return compute_inner_prod(A12_sigma, sigma);
|
|
|
+ };
|
|
|
+ auto compute_A22adj = [&compute_A22] (const Real alpha) {
|
|
|
+ return compute_A22(alpha);
|
|
|
+ };
|
|
|
+ auto compute_Aadj = [&compute_A11adj,&compute_A12adj,&compute_A21adj,&compute_A22adj] (const Vector<Real>& x) {
|
|
|
+ const Vector<Real> sigma(x.Dim()-1,(Iterator<Real>)x.begin(),false);
|
|
|
+ const Real& alpha = x[x.Dim()-1];
|
|
|
+
|
|
|
+ Vector<Real> Ax;
|
|
|
+ Ax.ReInit(x.Dim());
|
|
|
+ Vector<Real> Bdotn(x.Dim()-1,Ax.begin(),false);
|
|
|
+ Real& flux = Ax[x.Dim()-1];
|
|
|
+
|
|
|
+ Vector<Real> Adotn_0, Adotn_1;
|
|
|
+ compute_A11adj(Adotn_0, sigma);
|
|
|
+ compute_A21adj(Adotn_1, alpha);
|
|
|
+ Bdotn = Adotn_0 + Adotn_1;
|
|
|
+
|
|
|
+ flux = compute_A12adj(sigma) + compute_A22adj(alpha);
|
|
|
+ return Ax;
|
|
|
+ };
|
|
|
+ auto compute_invAadj = [&S,&comm,&compute_Aadj] (Vector<Real>& b) {
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&compute_Aadj](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
|
|
|
+ (*Ax) = compute_Aadj(x);
|
|
|
+ };
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<Real> x(b.Dim());
|
|
|
+ x = 0;
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&x, BIOp, b, 1e-8, 50);
|
|
|
+ return x;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_dg_dsigma = [&S, &normal, &compute_dot_prod](const Vector<ElemBasis>& B) { // dg_dsigma = \int 2 B \cdot (\nabla G + n/2)
|
|
|
+ Vector<ElemBasis> B_dot_gradG;
|
|
|
+ S.quadrature_dUxF.Eval(B_dot_gradG, S.GetElemList(), B, S.Laplace_dUxF);
|
|
|
+ return B_dot_gradG * (-2.0) + compute_dot_prod(B,normal);
|
|
|
+ };
|
|
|
+ auto compute_dg_dalpha = [&S,&compute_B0,&compute_inner_prod] (const Vector<ElemBasis>& B) {
|
|
|
+ auto dB_dalpha = compute_B0(1);
|
|
|
+ return 2*compute_inner_prod(B,dB_dalpha);
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_dg_dnu = [&S,&comm,&normal,&compute_inner_prod,&area_elem,&compute_dB0](const Vector<ElemBasis>& sigma, Real alpha, const Vector<ElemBasis>& B) { // dg_dnu = (B*B) 2H - (2 B) \cdot (n \cdnot nabla) \nabla G[sigma] + (2 B) \alpha dB0_dnu \hat{\theta} + sigma (\nabla D)^T [2 B] + (2H) sigma (\nabla G)^T [2 B]
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+ Vector<ElemBasis> v = B * 2.0;
|
|
|
+
|
|
|
+ Vector<ElemBasis> dg_dnu0(Nelem), dg_dnu1(Nelem), dg_dnu2(Nelem), dg_dnu3(Nelem), dg_dnu4(Nelem);
|
|
|
+ dg_dnu0 = 0;
|
|
|
+ dg_dnu1 = 0;
|
|
|
+ dg_dnu2 = 0;
|
|
|
+ dg_dnu3 = 0;
|
|
|
+ dg_dnu4 = 0;
|
|
|
+
|
|
|
+ Vector<ElemBasis> H(Nelem);
|
|
|
+ { // Set mean curvature H
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ Vector<ElemBasis> dX, d2X;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+ ElemBasis::Grad(d2X, dX);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,2,2> I, invI, II;
|
|
|
+ for (Long k0 = 0; k0 < 2; k0++) {
|
|
|
+ for (Long k1 = 0; k1 < 2; k1++) {
|
|
|
+ I(k0,k1) = 0;
|
|
|
+ I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j];
|
|
|
+ I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j];
|
|
|
+ I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j];
|
|
|
+
|
|
|
+ II(k0,k1) = 0;
|
|
|
+ II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ { // Set invI
|
|
|
+ Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0);
|
|
|
+ invI(0,0) = I(1,1) / detI;
|
|
|
+ invI(0,1) = -I(0,1) / detI;
|
|
|
+ invI(1,0) = -I(1,0) / detI;
|
|
|
+ invI(1,1) = I(0,0) / detI;
|
|
|
+ }
|
|
|
+ { // Set H
|
|
|
+ H[i][j] = 0;
|
|
|
+ H[i][j] += -0.5 * II(0,0)*invI(0,0);
|
|
|
+ H[i][j] += -0.5 * II(0,1)*invI(0,1);
|
|
|
+ H[i][j] += -0.5 * II(1,0)*invI(1,0);
|
|
|
+ H[i][j] += -0.5 * II(1,1)*invI(1,1);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dg_dnu = (B*B) 2H
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dnu0[i][j] = 0;
|
|
|
+ dg_dnu0[i][j] += B[i*COORD_DIM+0][j] * B[i*COORD_DIM+0][j] * (2.0*H[i][j]);
|
|
|
+ dg_dnu0[i][j] += B[i*COORD_DIM+1][j] * B[i*COORD_DIM+1][j] * (2.0*H[i][j]);
|
|
|
+ dg_dnu0[i][j] += B[i*COORD_DIM+2][j] * B[i*COORD_DIM+2][j] * (2.0*H[i][j]);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dg_dnu1 = (2 B) \cdot (n \cdnot nabla) \nabla G[sigma]
|
|
|
+ Vector<ElemBasis> d2Gsigma;
|
|
|
+ Quadrature<Real> quadrature_Fxd2U;
|
|
|
+ quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
|
|
|
+ quadrature_Fxd2U.Eval(d2Gsigma, S.GetElemList(), sigma, S.Laplace_Fxd2U);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dnu1[i][j] = 0;
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+0][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j];
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+1][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+1][j];
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+2][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+3][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+0][j];
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+4][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j];
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+5][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+6][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+0][j];
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+7][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+1][j];
|
|
|
+ dg_dnu1[i][j] -= d2Gsigma[i*9+8][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dg_dnu2 = (2 B) \alpha dB0_dnu \hat{\theta}
|
|
|
+ Vector<ElemBasis> dB0 = compute_dB0(alpha);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dnu2[i][j] = 0;
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+0][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j];
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+1][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+0][j];
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+2][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+0][j];
|
|
|
+
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+3][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+1][j];
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+4][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j];
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+5][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+1][j];
|
|
|
+
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+6][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+2][j];
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+7][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+2][j];
|
|
|
+ dg_dnu2[i][j] += dB0[i*9+8][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dg_dnu3 = (sigma (\nabla D)^T [2 B]
|
|
|
+ Vector<ElemBasis> nablaDtv;
|
|
|
+ Quadrature<Real> quadrature_dUxD;
|
|
|
+ quadrature_dUxD.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
|
|
|
+ quadrature_dUxD.Eval(nablaDtv, S.GetElemList(), v, S.Laplace_dUxD);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dnu3[i][j] = 0;
|
|
|
+ dg_dnu3[i][j] += sigma[i][j] * nablaDtv[i*COORD_DIM+0][j]*normal[i*COORD_DIM+0][j];
|
|
|
+ dg_dnu3[i][j] += sigma[i][j] * nablaDtv[i*COORD_DIM+1][j]*normal[i*COORD_DIM+1][j];
|
|
|
+ dg_dnu3[i][j] += sigma[i][j] * nablaDtv[i*COORD_DIM+2][j]*normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dg_dnu4 = (2H) sigma (\nabla G)^T [2 B]
|
|
|
+ Quadrature<Real> quadrature_dUxF;
|
|
|
+ quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
|
|
|
+ quadrature_dUxF.Eval(dg_dnu4, S.GetElemList(), v, S.Laplace_dUxF);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dnu4[i][j] *= 2*H[i][j] * sigma[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ return dg_dnu0 + dg_dnu1 + dg_dnu2 + dg_dnu3 - dg_dnu4;
|
|
|
+ };
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ Real flux = 1.0, alpha;
|
|
|
+ Vector<ElemBasis> sigma;
|
|
|
+ compute_invA(sigma, alpha, flux);
|
|
|
+ Vector<ElemBasis> B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
|
|
|
+ Real g = compute_inner_prod(B, B);
|
|
|
+ std::cout<<"g = "<<g<<'\n';
|
|
|
+
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), sigma, ORDER);
|
|
|
+ vtu.WriteVTK("sigma", comm);
|
|
|
+ }
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), B, ORDER);
|
|
|
+ vtu.WriteVTK("B", comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ if (0) { // test dg_dnu
|
|
|
+ auto compute_g = [&S,&comm,&normal,&area_elem,&sigma,&alpha,&compute_norm_area_elem,&compute_B0,&compute_inner_prod](const Vector<ElemBasis>& nu, Real eps) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> X_orig(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ X_orig[i*COORD_DIM+0][j] = S.Elem(i,0)[j];
|
|
|
+ X_orig[i*COORD_DIM+1][j] = S.Elem(i,1)[j];
|
|
|
+ X_orig[i*COORD_DIM+2][j] = S.Elem(i,2)[j];
|
|
|
+ S.Elem(i,0)[j] += eps*nu[i][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] += eps*nu[i][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+
|
|
|
+ Vector<Real> Xt(Nelem*Nnodes*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ for (Long k = 0; k < COORD_DIM; k++) {
|
|
|
+ Xt[(i*Nnodes+j)*COORD_DIM+k] = S.Elem(i,k)[j] - 1e-4*normal[i*COORD_DIM+k][j];// + eps*nu[i][j] * normal[i*COORD_DIM+k][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Vector<ElemBasis> B0 = compute_B0(alpha);
|
|
|
+
|
|
|
+ Vector<ElemBasis> B1;
|
|
|
+ Quadrature<Real> quadrature_FxdU;
|
|
|
+ quadrature_FxdU.template Setup<ElemBasis>(S.GetElemList(), Xt, S.Laplace_FxdU, order_singular, order_direct, -1, comm);
|
|
|
+ quadrature_FxdU.Eval(B1, S.GetElemList(), sigma, S.Laplace_FxdU);
|
|
|
+
|
|
|
+ Real g = compute_inner_prod(B0+B1, B0+B1);
|
|
|
+
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ S.Elem(i,0)[j] = X_orig[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] = X_orig[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+
|
|
|
+ return g;
|
|
|
+ };
|
|
|
+
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> nu(Nelem);
|
|
|
+ nu = 1; //area_elem;
|
|
|
+
|
|
|
+ Vector<ElemBasis> dg_dnu = compute_dg_dnu(sigma, alpha, B);
|
|
|
+ std::cout<<compute_inner_prod(dg_dnu, nu)<<'\n';
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), dg_dnu, ORDER);
|
|
|
+ vtu.WriteVTK("dg_dnu", comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ Real eps = 1e-5;
|
|
|
+ Real g0 = compute_g(nu,-eps);
|
|
|
+ Real g1 = compute_g(nu,eps);
|
|
|
+ std::cout<<"g = "<<g0<<" g = "<<g1<<" dg_dnu = "<<(g1-g0)/(2*eps)<<'\n';
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ if (0) { // test dg_dsigma
|
|
|
+ Vector<ElemBasis> dg_dsigma = compute_dg_dsigma(B);
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), dg_dsigma, ORDER);
|
|
|
+ vtu.WriteVTK("dg_dsigma", comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ Real dt = 1e-1;
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const auto& quad_wts = ElemBasis::QuadWts();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dg_dsigma_(Nelem);
|
|
|
+ dg_dsigma_ = 0;
|
|
|
+ for (Long i = 0; i < Nelem; i++) { // Set dg_dsigma_
|
|
|
+ for (Long j = 0; j < ElemBasis::Size(); j++) {
|
|
|
+ auto sigma_0 = sigma;
|
|
|
+ auto sigma_1 = sigma;
|
|
|
+ sigma_0[i][j] -= 0.5*dt;
|
|
|
+ sigma_1[i][j] += 0.5*dt;
|
|
|
+ auto B_0 = compute_half_n_plus_dG(sigma_0) + compute_B0(alpha);
|
|
|
+ auto B_1 = compute_half_n_plus_dG(sigma_1) + compute_B0(alpha);
|
|
|
+ auto g_0 = compute_inner_prod(B_0, B_0);
|
|
|
+ auto g_1 = compute_inner_prod(B_1, B_1);
|
|
|
+ dg_dsigma_[i][j] = (g_1 - g_0) / dt;
|
|
|
+ dg_dsigma_[i][j] /= quad_wts[j] * area_elem[i][j];
|
|
|
+ std::cout<<dg_dsigma_[i][j]<<' '<<j<<' '<<ElemBasis::Size()<<'\n'; ////////////////
|
|
|
+ }
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), dg_dsigma_, ORDER);
|
|
|
+ vtu.WriteVTK("dg_dsigma_", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (0) { // test dg_dalpha
|
|
|
+ Real dg_dalpha = compute_dg_dalpha(B);
|
|
|
+
|
|
|
+ Real dt = 1e-1;
|
|
|
+ auto B_0 = compute_half_n_plus_dG(sigma) + compute_B0(alpha - 0.5*dt);
|
|
|
+ auto B_1 = compute_half_n_plus_dG(sigma) + compute_B0(alpha + 0.5*dt);
|
|
|
+ auto g_0 = compute_inner_prod(B_0, B_0);
|
|
|
+ auto g_1 = compute_inner_prod(B_1, B_1);
|
|
|
+ Real dg_dalpha_ = (g_1 - g_0) / dt;
|
|
|
+ std::cout<<dg_dalpha<<' '<<dg_dalpha_<<'\n';
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ if (0) { // test compute_A21adj
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<Real> A21adj_;
|
|
|
+ compute_A21adj(A21adj_, flux);
|
|
|
+ Vector<ElemBasis> A21adj(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ A21adj[i][j] = A21adj_[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), A21adj, ORDER);
|
|
|
+ vtu.WriteVTK("A21adj", comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ { // verify
|
|
|
+ Vector<Real> sigma_(Nelem*Nnodes);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i*Nnodes+j] = sigma[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Real flux = compute_inner_prod(A21adj, sigma);
|
|
|
+ std::cout<<"Error: "<<compute_A21(sigma_)-flux<<'\n';
|
|
|
+ }
|
|
|
+ { // compute finite-difference matrix
|
|
|
+ Real dt = 1e+1;
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const auto& quad_wts = ElemBasis::QuadWts();
|
|
|
+
|
|
|
+ Vector<ElemBasis> A21(Nelem);
|
|
|
+ A21 = 0;
|
|
|
+ for (Long i = 0; i < Nelem; i++) { // Set A21
|
|
|
+ for (Long j = 0; j < ElemBasis::Size(); j++) {
|
|
|
+ Vector<Real> sigma_0(Nelem*ElemBasis::Size());
|
|
|
+ Vector<Real> sigma_1(Nelem*ElemBasis::Size());
|
|
|
+ sigma_0 = 0;
|
|
|
+ sigma_1 = 0;
|
|
|
+ sigma_0[i*ElemBasis::Size()+j] -= 0.5*dt;
|
|
|
+ sigma_1[i*ElemBasis::Size()+j] += 0.5*dt;
|
|
|
+ auto flux_0 = compute_A21(sigma_0);
|
|
|
+ auto flux_1 = compute_A21(sigma_1);
|
|
|
+ A21[i][j] = (flux_1 - flux_0) / dt;
|
|
|
+ A21[i][j] /= quad_wts[j] * area_elem[i][j];
|
|
|
+ std::cout<<A21[i][j]<<' '<<j<<' '<<ElemBasis::Size()<<'\n'; ////////////////
|
|
|
+ }
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), A21, ORDER);
|
|
|
+ vtu.WriteVTK("A21", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ auto compute_invA11 = [&S,&normal,&comm,&compute_A11](const Vector<ElemBasis>& rhs) { // Solver for sigma: sigma/2 + n.dG[sigma] = rhs
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&S,&normal,&compute_A11](sctl::Vector<Real>* A11_sigma, const sctl::Vector<Real>& sigma) {
|
|
|
+ compute_A11(*A11_sigma, sigma);
|
|
|
+ };
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> sigma(Nelem);
|
|
|
+ Vector<Real> rhs_(Nelem * Nnodes), sigma_(Nelem * Nnodes);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ rhs_[i*Nnodes+j] = rhs[i][j];
|
|
|
+ sigma_[i*Nnodes+j] = 0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&sigma_, BIOp, rhs_, 1e-8, 50);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma[i][j] = sigma_[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return sigma;
|
|
|
+ };
|
|
|
+ auto compute_invA11adj = [&S,&normal,&comm,&compute_A11adj](const Vector<ElemBasis>& rhs) { // Solver for sigma: A11adj sigma = rhs
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&S,&compute_A11adj](sctl::Vector<Real>* A11adj_sigma, const sctl::Vector<Real>& sigma) {
|
|
|
+ compute_A11adj(*A11adj_sigma, sigma);
|
|
|
+ };
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> sigma(Nelem);
|
|
|
+ Vector<Real> rhs_(Nelem * Nnodes), sigma_(Nelem * Nnodes);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ rhs_[i*Nnodes+j] = rhs[i][j];
|
|
|
+ sigma_[i*Nnodes+j] = 0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&sigma_, BIOp, rhs_, 1e-8, 50);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma[i][j] = sigma_[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return sigma;
|
|
|
+ };
|
|
|
+ if (0) { // Test invA11adj
|
|
|
+ Vector<ElemBasis> dg_dsigma = compute_dg_dsigma(B);
|
|
|
+ Real a = compute_inner_prod(dg_dsigma, compute_invA11(sigma));
|
|
|
+ Real b = compute_inner_prod(compute_invA11adj(dg_dsigma), sigma);
|
|
|
+ std::cout<<a<<' '<<b<<'\n';
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ // 0.168275 0.117983 -0.110446 -96.7293
|
|
|
+ // 0.603869 -1.901900 -1.229930 -245.5050
|
|
|
+ auto compute_u_dAdnu_v_00 = [&S,&normal,&comm,&compute_half_n_plus_dG,&compute_grad_adj] (const Vector<Real>& u_, const Vector<Real>& v_) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> u(Nelem), u_n(Nelem*COORD_DIM), v(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ u[i][j] = u_[i*Nnodes+j];
|
|
|
+ v[i][j] = v_[i*Nnodes+j];
|
|
|
+ u_n[i*COORD_DIM+0][j] = u[i][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ u_n[i*COORD_DIM+1][j] = u[i][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ u_n[i*COORD_DIM+2][j] = u[i][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> dAdnu0(Nelem), dAdnu1(Nelem), dAdnu2(Nelem), dAdnu3(Nelem);
|
|
|
+ dAdnu0 = 0;
|
|
|
+ dAdnu1 = 0;
|
|
|
+ dAdnu2 = 0;
|
|
|
+ dAdnu3 = 0;
|
|
|
+
|
|
|
+ Vector<ElemBasis> H(Nelem);
|
|
|
+ { // Set mean curvature H
|
|
|
+ const Vector<ElemBasis> X = S.GetElemList().ElemVector();
|
|
|
+ Vector<ElemBasis> dX, d2X;
|
|
|
+ ElemBasis::Grad(dX, X);
|
|
|
+ ElemBasis::Grad(d2X, dX);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Tensor<Real,true,2,2> I, invI, II;
|
|
|
+ for (Long k0 = 0; k0 < 2; k0++) {
|
|
|
+ for (Long k1 = 0; k1 < 2; k1++) {
|
|
|
+ I(k0,k1) = 0;
|
|
|
+ I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j];
|
|
|
+ I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j];
|
|
|
+ I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j];
|
|
|
+
|
|
|
+ II(k0,k1) = 0;
|
|
|
+ II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ { // Set invI
|
|
|
+ Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0);
|
|
|
+ invI(0,0) = I(1,1) / detI;
|
|
|
+ invI(0,1) = -I(0,1) / detI;
|
|
|
+ invI(1,0) = -I(1,0) / detI;
|
|
|
+ invI(1,1) = I(0,0) / detI;
|
|
|
+ }
|
|
|
+ { // Set H
|
|
|
+ H[i][j] = 0;
|
|
|
+ H[i][j] += -0.5 * II(0,0)*invI(0,0);
|
|
|
+ H[i][j] += -0.5 * II(0,1)*invI(0,1);
|
|
|
+ H[i][j] += -0.5 * II(1,0)*invI(1,0);
|
|
|
+ H[i][j] += -0.5 * II(1,1)*invI(1,1);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dAdnu0 = u B \cdot grad_nu
|
|
|
+ Vector<ElemBasis> B = compute_half_n_plus_dG(v);
|
|
|
+ Vector<ElemBasis> u_B(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ u_B[i*COORD_DIM+0][j] = u[i][j] * B[i*COORD_DIM+0][j];
|
|
|
+ u_B[i*COORD_DIM+1][j] = u[i][j] * B[i*COORD_DIM+1][j];
|
|
|
+ u_B[i*COORD_DIM+2][j] = u[i][j] * B[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ dAdnu0 = compute_grad_adj(u_B);
|
|
|
+
|
|
|
+ // dAdnu1 = (2H) v (I/2 + \nabla G)^T [u n]
|
|
|
+ Quadrature<Real> quadrature_dUxF;
|
|
|
+ quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
|
|
|
+ quadrature_dUxF.Eval(dAdnu1, S.GetElemList(), u_n, S.Laplace_dUxF);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dAdnu1[i][j] *= -2*H[i][j] * v[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dAdnu2 = (u n) \cdot (n \cdnot \nabla) \nabla G[v]
|
|
|
+ Vector<ElemBasis> d2G_v;
|
|
|
+ Quadrature<Real> quadrature_Fxd2U;
|
|
|
+ quadrature_Fxd2U.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_Fxd2U, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
|
|
|
+ quadrature_Fxd2U.Eval(d2G_v, S.GetElemList(), v, S.Laplace_Fxd2U);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dAdnu2[i][j] = 0;
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+0][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+0][j];
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+1][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+1][j];
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+2][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+3][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+0][j];
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+4][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+1][j];
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+5][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+6][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+0][j];
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+7][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+1][j];
|
|
|
+ dAdnu2[i][j] -= d2G_v[i*9+8][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ // dAdnu3 = (v (\nabla D)^T [u n]
|
|
|
+ Vector<ElemBasis> nablaDt_u_n;
|
|
|
+ Quadrature<Real> quadrature_dUxD;
|
|
|
+ quadrature_dUxD.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxD, order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER));
|
|
|
+ quadrature_dUxD.Eval(nablaDt_u_n, S.GetElemList(), u_n, S.Laplace_dUxD);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dAdnu3[i][j] = 0;
|
|
|
+ dAdnu3[i][j] += v[i][j] * nablaDt_u_n[i*COORD_DIM+0][j]*normal[i*COORD_DIM+0][j];
|
|
|
+ dAdnu3[i][j] += v[i][j] * nablaDt_u_n[i*COORD_DIM+1][j]*normal[i*COORD_DIM+1][j];
|
|
|
+ dAdnu3[i][j] += v[i][j] * nablaDt_u_n[i*COORD_DIM+2][j]*normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ return dAdnu0 + dAdnu1 + dAdnu2 + dAdnu3;
|
|
|
+ };
|
|
|
+ auto compute_u_dAdnu_v_01 = [&S,&comm,&compute_dB0,&normal,&area_elem,&compute_B0,&compute_grad_adj] (const Vector<Real>& u, const Vector<Real>& v) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dAdnu(Nelem);
|
|
|
+ Vector<ElemBasis> dB0 = compute_dB0(v[Nelem*Nnodes]);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Real n_n_dB0 = 0;
|
|
|
+ n_n_dB0 += dB0[i*9+0][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ n_n_dB0 += dB0[i*9+1][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ n_n_dB0 += dB0[i*9+2][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j];
|
|
|
+
|
|
|
+ n_n_dB0 += dB0[i*9+3][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ n_n_dB0 += dB0[i*9+4][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ n_n_dB0 += dB0[i*9+5][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j];
|
|
|
+
|
|
|
+ n_n_dB0 += dB0[i*9+6][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ n_n_dB0 += dB0[i*9+7][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ n_n_dB0 += dB0[i*9+8][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
|
|
|
+
|
|
|
+ dAdnu[i][j] = u[i*Nnodes+j] * n_n_dB0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> B0 = compute_B0(v[Nelem*Nnodes]);
|
|
|
+ Vector<ElemBasis> u_B0(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ u_B0[i*COORD_DIM+0][j] = u[i*Nnodes+j] * B0[i*COORD_DIM+0][j];
|
|
|
+ u_B0[i*COORD_DIM+1][j] = u[i*Nnodes+j] * B0[i*COORD_DIM+1][j];
|
|
|
+ u_B0[i*COORD_DIM+2][j] = u[i*Nnodes+j] * B0[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ dAdnu += compute_grad_adj(u_B0);
|
|
|
+
|
|
|
+ return dAdnu;
|
|
|
+ };
|
|
|
+ auto compute_u_dAdnu_v_10 = [&S,&comm] (const Vector<Real>& u, const Vector<Real>& v) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dAdnu(Nelem*Nnodes);
|
|
|
+ dAdnu = 0;
|
|
|
+ return dAdnu;
|
|
|
+ };
|
|
|
+ auto compute_u_dAdnu_v_11 = [&S,&comm] (const Vector<Real>& u, const Vector<Real>& v) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> dAdnu(Nelem*Nnodes);
|
|
|
+ dAdnu = 0;
|
|
|
+ return dAdnu;
|
|
|
+ };
|
|
|
+
|
|
|
+ auto compute_Av = [&S,&area_elem,&normal,&compute_norm_area_elem,&compute_A,&comm] (const Vector<Real>& v, const Vector<ElemBasis>& nu, Real eps) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> X_orig(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ X_orig[i*COORD_DIM+0][j] = S.Elem(i,0)[j];
|
|
|
+ X_orig[i*COORD_DIM+1][j] = S.Elem(i,1)[j];
|
|
|
+ X_orig[i*COORD_DIM+2][j] = S.Elem(i,2)[j];
|
|
|
+ S.Elem(i,0)[j] += eps*nu[i][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] += eps*nu[i][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ Vector<Real> Av = compute_A(v);
|
|
|
+
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ S.Elem(i,0)[j] = X_orig[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] = X_orig[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ return Av;
|
|
|
+ };
|
|
|
+ auto compute_u_dAdnu_v = [&S,&compute_Av,&compute_inner_prod] (const Vector<Real>& u, const Vector<Real>& v, const Vector<ElemBasis>& nu) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Real eps = 1e-5;
|
|
|
+ Vector<Real> Av0 = compute_Av(v,nu,-eps);
|
|
|
+ Vector<Real> Av1 = compute_Av(v,nu,eps);
|
|
|
+ Vector<Real> dAdnu_v = (Av1-Av0)*(1/(2*eps));
|
|
|
+
|
|
|
+ Real u_dAdnu_v;
|
|
|
+ { // set u_dAdnu_v
|
|
|
+ Vector<ElemBasis> u_(Nelem), dAdnu_v_(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ u_[i][j] = u[i*Nnodes+j];
|
|
|
+ dAdnu_v_[i][j] = dAdnu_v[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ u_dAdnu_v = compute_inner_prod(u_, dAdnu_v_);
|
|
|
+ u_dAdnu_v += u[Nelem*Nnodes] * dAdnu_v[Nelem*Nnodes];
|
|
|
+ }
|
|
|
+ return u_dAdnu_v;
|
|
|
+ };
|
|
|
+ if (0) { // test dA_dnu
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> nu(Nelem);
|
|
|
+ Vector<Real> u(Nelem*Nnodes+1), v(Nelem*Nnodes+1);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ v[i*Nnodes+j] = sigma[i][j];
|
|
|
+ u[i*Nnodes+j] = sigma[i][j]*area_elem[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ v[Nelem*Nnodes] = 0; //alpha;
|
|
|
+ u[Nelem*Nnodes] = 0;
|
|
|
+ nu = 1; //area_elem;
|
|
|
+
|
|
|
+ Real u_dAdnu_v = compute_u_dAdnu_v(u, v, nu);
|
|
|
+ std::cout<<"u_dAdnu_v = "<<u_dAdnu_v<<'\n';
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ auto compute_dsigma_dnu = [&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&comm] (const Vector<ElemBasis>& nu, Real eps) {
|
|
|
+ auto compute_sigma = [&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&comm] (const Vector<ElemBasis>& nu, Real eps) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> X_orig(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ X_orig[i*COORD_DIM+0][j] = S.Elem(i,0)[j];
|
|
|
+ X_orig[i*COORD_DIM+1][j] = S.Elem(i,1)[j];
|
|
|
+ X_orig[i*COORD_DIM+2][j] = S.Elem(i,2)[j];
|
|
|
+ S.Elem(i,0)[j] += eps*nu[i][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] += eps*nu[i][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ Real flux = 1.0, alpha;
|
|
|
+ Vector<ElemBasis> sigma;
|
|
|
+ compute_invA(sigma, alpha, flux);
|
|
|
+ Vector<Real> sigma_(Nelem*Nnodes+1);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i*Nnodes+j] = sigma[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ sigma_[Nelem*Nnodes] = alpha;
|
|
|
+
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ S.Elem(i,0)[j] = X_orig[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] = X_orig[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ return sigma_;
|
|
|
+ };
|
|
|
+ auto sigma0 = compute_sigma(nu,-eps);
|
|
|
+ auto sigma1 = compute_sigma(nu,eps);
|
|
|
+ return (sigma1-sigma0) * (1/(2*eps));
|
|
|
+ };
|
|
|
+ if (0) { // verify dA_dnu sigma + A dsigma_dnu = 0
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> nu(Nelem);
|
|
|
+ nu = 1; //area_elem;
|
|
|
+
|
|
|
+ Vector<Real> dA_dnu_sigma;
|
|
|
+ { // Set dA_dnu_simga
|
|
|
+ Vector<Real> sigma_(Nelem*Nnodes+1);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i*Nnodes+j] = sigma[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ sigma_[Nelem*Nnodes] = alpha;
|
|
|
+
|
|
|
+ Real eps = 1e-3;
|
|
|
+ Vector<Real> Asigma0 = compute_Av(sigma_,nu,-eps);
|
|
|
+ Vector<Real> Asigma1 = compute_Av(sigma_,nu,eps);
|
|
|
+ dA_dnu_sigma = (Asigma1-Asigma0) * (1/(2*eps));
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<Real> A_dsigma_dnu;
|
|
|
+ { // Set A_dsigma_dnu
|
|
|
+ Vector<Real> dsigma_dnu = compute_dsigma_dnu(nu, 1e-3);
|
|
|
+ A_dsigma_dnu = compute_A(dsigma_dnu);
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<ElemBasis> dA_dnu_sigma_(Nelem);
|
|
|
+ Vector<ElemBasis> A_dsigma_dnu_(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dA_dnu_sigma_[i][j] = dA_dnu_sigma[i*Nnodes+j];
|
|
|
+ A_dsigma_dnu_[i][j] = A_dsigma_dnu[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ std::cout<<dA_dnu_sigma[Nelem*Nnodes] + A_dsigma_dnu[Nelem*Nnodes]<<'\n';
|
|
|
+
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu;
|
|
|
+ vtu.AddElems(S.GetElemList(), dA_dnu_sigma_ + A_dsigma_dnu_, ORDER);
|
|
|
+ vtu.WriteVTK("err", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ if (1) { // test grad_g
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> nu(Nelem);
|
|
|
+ nu = 1; //area_elem;
|
|
|
+
|
|
|
+ if (1) {
|
|
|
+ Real dg_dnu0 = compute_inner_prod(nu, compute_dg_dnu(sigma, alpha, B));
|
|
|
+
|
|
|
+ Vector<Real> dg_dsigma(Nelem*Nnodes+1);
|
|
|
+ { // Set dg_dsigma
|
|
|
+ Vector<ElemBasis> dg_dsigma_ = compute_dg_dsigma(B);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dsigma[i*Nnodes+j] = dg_dsigma_[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ dg_dsigma[Nelem*Nnodes] = compute_dg_dalpha(B);
|
|
|
+ }
|
|
|
+
|
|
|
+ Real dg_dnu1, dg_dnu2, dg_dnu3, dg_dnu4;
|
|
|
+ if (1) { // Set dg_dnu = - dg_dsigma invA dA_dnu sigma
|
|
|
+ Vector<Real> dg_dsigma_invA = compute_invAadj(dg_dsigma);
|
|
|
+ Vector<Real> sigma_(Nelem*Nnodes+1);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i*Nnodes+j] = sigma[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ sigma_[Nelem*Nnodes] = alpha;
|
|
|
+
|
|
|
+ dg_dnu1 = -compute_inner_prod(nu, compute_u_dAdnu_v_00(dg_dsigma_invA, sigma_));
|
|
|
+ dg_dnu2 = -compute_inner_prod(nu, compute_u_dAdnu_v_01(dg_dsigma_invA, sigma_));
|
|
|
+ dg_dnu3 = -compute_inner_prod(nu, compute_u_dAdnu_v_10(dg_dsigma_invA, sigma_));
|
|
|
+ dg_dnu4 = -compute_inner_prod(nu, compute_u_dAdnu_v_11(dg_dsigma_invA, sigma_));
|
|
|
+ std::cout<<dg_dnu1<<' '<<dg_dnu2<<' '<<dg_dnu3<<' '<<dg_dnu4<<'\n';
|
|
|
+ exit(0);
|
|
|
+ }
|
|
|
+ if (0) { // Set dg_dnu = - dg_dsigma invA dA_dnu sigma
|
|
|
+ Vector<Real> dg_dsigma_invA = compute_invAadj(dg_dsigma);
|
|
|
+ Vector<Real> sigma_(Nelem*Nnodes+1);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ sigma_[i*Nnodes+j] = sigma[i][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ sigma_[Nelem*Nnodes] = alpha;
|
|
|
+
|
|
|
+ Vector<Real> dg_dsigma_invA_0 = dg_dsigma_invA; dg_dsigma_invA_0[Nelem*Nnodes] = 0;
|
|
|
+ Vector<Real> dg_dsigma_invA_1(Nelem*Nnodes+1); dg_dsigma_invA_1 = 0; dg_dsigma_invA_1[Nelem*Nnodes] = dg_dsigma_invA[Nelem*Nnodes];
|
|
|
+
|
|
|
+ Vector<Real> sigma_0 = sigma_; sigma_0[Nelem*Nnodes] = 0;
|
|
|
+ Vector<Real> sigma_1(Nelem*Nnodes+1); sigma_1 = 0; sigma_1[Nelem*Nnodes] = sigma_[Nelem*Nnodes];
|
|
|
+
|
|
|
+ dg_dnu1 = -compute_u_dAdnu_v(dg_dsigma_invA_0, sigma_0, nu);
|
|
|
+ dg_dnu2 = -compute_u_dAdnu_v(dg_dsigma_invA_0, sigma_1, nu);
|
|
|
+ dg_dnu3 = -compute_u_dAdnu_v(dg_dsigma_invA_1, sigma_0, nu);
|
|
|
+ dg_dnu4 = -compute_u_dAdnu_v(dg_dsigma_invA_1, sigma_1, nu);
|
|
|
+ std::cout<<dg_dnu1<<' '<<dg_dnu2<<' '<<dg_dnu3<<' '<<dg_dnu4<<'\n';
|
|
|
+ }
|
|
|
+ if (0) { // Set dg_dnu = dg_dsigma dsigma_dnu
|
|
|
+ Vector<Real> dsigma_dnu = compute_dsigma_dnu(nu, 1e-3);
|
|
|
+ Vector<ElemBasis> dg_dsigma_(Nelem), dsigma_dnu_(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ dg_dsigma_[i][j] = dg_dsigma[i*Nnodes+j];
|
|
|
+ dsigma_dnu_[i][j] = dsigma_dnu[i*Nnodes+j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ dg_dnu1 = compute_inner_prod(dg_dsigma_, dsigma_dnu_);
|
|
|
+ dg_dnu1 += dg_dsigma[Nelem*Nnodes] * dsigma_dnu[Nelem*Nnodes];
|
|
|
+ dg_dnu2 = 0;
|
|
|
+ dg_dnu3 = 0;
|
|
|
+ dg_dnu4 = 0;
|
|
|
+ }
|
|
|
+
|
|
|
+ std::cout<<dg_dnu0<<' '<<dg_dnu1+dg_dnu2+dg_dnu3+dg_dnu4<<'\n';
|
|
|
+ }
|
|
|
+
|
|
|
+ auto compute_g = [&sigma,&alpha,&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&compute_half_n_plus_dG,&compute_B0,&compute_inner_prod,&comm] (const Vector<ElemBasis>& nu, Real eps) {
|
|
|
+ const Long Nelem = S.GetElemList().NElem();
|
|
|
+ const Long Nnodes = ElemBasis::Size();
|
|
|
+
|
|
|
+ Vector<ElemBasis> X_orig(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ X_orig[i*COORD_DIM+0][j] = S.Elem(i,0)[j];
|
|
|
+ X_orig[i*COORD_DIM+1][j] = S.Elem(i,1)[j];
|
|
|
+ X_orig[i*COORD_DIM+2][j] = S.Elem(i,2)[j];
|
|
|
+ S.Elem(i,0)[j] += eps*nu[i][j] * normal[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] += eps*nu[i][j] * normal[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ Real flux = 1.0, alpha;
|
|
|
+ Vector<ElemBasis> sigma;
|
|
|
+ compute_invA(sigma, alpha, flux);
|
|
|
+ Vector<ElemBasis> B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
|
|
|
+ Real g = compute_inner_prod(B, B);
|
|
|
+
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ S.Elem(i,0)[j] = X_orig[i*COORD_DIM+0][j];
|
|
|
+ S.Elem(i,1)[j] = X_orig[i*COORD_DIM+1][j];
|
|
|
+ S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ compute_norm_area_elem(normal, area_elem);
|
|
|
+ S.quadrature_FxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU .template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_DxU , order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxdU.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_FxdU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_dUxF.template Setup<ElemBasis, ElemBasis>(S.GetElemList(), S.Laplace_dUxF, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ return g;
|
|
|
+ };
|
|
|
+ Real eps = 1e-3;
|
|
|
+ Real g0 = compute_g(nu,-eps);
|
|
|
+ Real g1 = compute_g(nu,eps);
|
|
|
+ std::cout<<"g = "<<g0<<" g = "<<g1<<" dg_dnu = "<<(g1-g0)/(2*eps)<<'\n';
|
|
|
+ }
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ // dg_dnu
|
|
|
+
|
|
|
+ // dA_dnu_sigma
|
|
|
+
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
+
|
|
|
+ //Profile::print(&comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ private:
|
|
|
+ void InitSurf(Long l) {
|
|
|
+ const auto& nodes = ElemBasis::Nodes();
|
|
|
+ const Long Nt = NtNp_[l*2+0];
|
|
|
+ const Long Np = NtNp_[l*2+1];
|
|
|
+
|
|
|
+ for (Long i = 0; i < Nt; i++) {
|
|
|
+ for (Long j = 0; j < Np; j++) {
|
|
|
+ for (Long k = 0; k < ElemBasis::Size(); k++) {
|
|
|
+ Real theta = (i + nodes[0][k]) * 2*const_pi<Real>()/Nt;
|
|
|
+ Real phi = (j + nodes[1][k]) * 2*const_pi<Real>()/Np;
|
|
|
+ Real X,Y,Z;
|
|
|
+ SurfGeom(X,Y,Z,theta,phi);
|
|
|
+ Elem(ElemIdx(l,i,j),0)[k] = X;
|
|
|
+ Elem(ElemIdx(l,i,j),1)[k] = Y;
|
|
|
+ Elem(ElemIdx(l,i,j),2)[k] = Z;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ static void SurfGeom(Real& X, Real& Y, Real& Z, Real theta, Real phi) {
|
|
|
+ sctl::Integer Nperiod = 5;
|
|
|
+#if 0
|
|
|
+ Real Aspect_ratio = 10.27932548522949;
|
|
|
+ Real coeffmat[21][21] = { 0.00000478813217, 0.00000000000000, 0.00000351611652, 0.00000135354389, 0.00000061357832, 0.00000220091101, 0.00000423862912, -0.00003000058678, 0.00000064187111, -0.00024228452821, 0.00003116775770, 0.00000176210710, 0.00000289141326, -0.00000150300525, 0.00000772853855, 0.00000098855242, 0.00000316606793, 0.00000002168364, 0.00000212047939, 0.00000299016097, 0.00000443224508,
|
|
|
+ 0.00000028202930, 0.00000000000000, -0.00000249222421, -0.00000203136278, 0.00000131104809, 0.00000011987446, -0.00000370760154, 0.00004553918916, -0.00007711342914, -0.00004685295062, 0.00011049838213, -0.00000197486270, 0.00000395827146, 0.00000615046474, 0.00000755337123, 0.00000700606006, 0.00000922725030, -0.00000043310337, 0.00000107416383, 0.00000449787694, 0.00000305137178,
|
|
|
+ 0.00001226376662, 0.00000000000000, 0.00000270820692, 0.00000208059305, 0.00000521478523, 0.00001779037302, 0.00000846544117, 0.00001120913385, -0.00065816845745, -0.00085107452469, -0.00013171190221, -0.00005540943675, -0.00001835885450, 0.00000101879823, 0.00000209222071, 0.00000091532502, -0.00000521515358, -0.00000209227142, -0.00000678545939, -0.00000034963549, -0.00000015111488,
|
|
|
+ 0.00001560274177, 0.00000000000000, 0.00000350691471, -0.00001160475040, -0.00001763036562, 0.00003487367940, -0.00002787247831, -0.00000910982726, 0.00008818832430, -0.00524408789352, 0.00009378376126, 0.00004184526188, 0.00002849263365, -0.00002757280527, 0.00003388467667, 0.00000706207265, 0.00000625263419, -0.00003315929280, -0.00001181772132, 0.00000311426015, 0.00001875682574,
|
|
|
+ -0.00000398287420, 0.00000000000000, -0.00001524541040, 0.00001724056165, 0.00002245173346, 0.00002806861812, -0.00000388776925, 0.00008143573359, -0.00005900909309, 0.00110496615525, 0.00134626252111, 0.00005128383054, -0.00001372421866, 0.00003612563887, 0.00002236580076, -0.00002728391883, 0.00001981237256, 0.00000655450458, 0.00000985319002, 0.00001347597299, 0.00000645987802,
|
|
|
+ 0.00003304968050, 0.00000000000000, -0.00000530822217, 0.00001324870937, -0.00003610889689, -0.00005478735329, -0.00005818806312, -0.00037112057908, -0.00017812002625, -0.00093204283621, 0.00115969858598, -0.00033559172880, -0.00010441876657, -0.00001617923044, -0.00000555065844, 0.00007343527250, -0.00004408047607, 0.00000403802142, 0.00001843931204, 0.00001694047933, 0.00001213414362,
|
|
|
+ -0.00000751115658, 0.00000000000000, 0.00005457974839, -0.00000334614515, 0.00005845565465, 0.00015000770509, 0.00021849104087, 0.00002724147635, 0.00167233624961, 0.00011666602222, 0.00276563479565, -0.00085952825611, -0.00030217235326, -0.00008841593808, 0.00000997664119, -0.00015285826521, 0.00002517224675, 0.00003009161810, 0.00001883217556, 0.00002146127554, 0.00001822445302,
|
|
|
+ -0.00004128706860, 0.00000000000000, -0.00003496417776, 0.00001088761655, -0.00000298955979, -0.00005359326315, -0.00019021633489, -0.00017992728681, -0.00347794801928, 0.00064632791327, 0.00449698418379, -0.00017710507382, 0.00006126180233, 0.00018059254216, 0.00002354096432, 0.00008189838991, -0.00010060678323, -0.00017183290038, 0.00019413756672, 0.00021334811754, 0.00011263617489,
|
|
|
+ 0.00000853522670, -0.00000000000000, -0.00006544789358, 0.00005424076880, -0.00000679056529, -0.00001249735487, -0.00053082982777, 0.00035396864405, -0.00115020677913, 0.05894451215863, 0.06573092192411, 0.01498018857092, 0.00278125284240, 0.00145188067108, 0.00033717858605, 0.00000800427370, -0.00009335305367, 0.00024286781263, -0.00023916347709, 0.00031213948387, 0.00018134393031,
|
|
|
+ -0.00002521496390, -0.00000000000000, -0.00054337945767, 0.00012690725271, 0.00053313979879, 0.00064233405283, -0.00047686311882, 0.00176536326762, 0.00074157933705, -0.02684566564858, 1.00000000000000, 0.07176169008017, 0.00837037432939, -0.00000381640211, 0.00088998704450, -0.00049218931235, -0.00024546548957, -0.00036608282244, 0.00049480766756, 0.00031158892671, 0.00006898906577,
|
|
|
+ 0.00021280418150, 0.00028127161204, -0.00070030166535, 0.00022237010126, -0.00028713891516, -0.00013800295710, 0.00005912094275, 0.00172126013786, -0.00618684850633, 0.03608432412148, Aspect_ratio , 0.49896776676178, 0.00091372377938, -0.00085712829605, -0.00124801427592, -0.00007427225501, -0.00005245858847, 0.00002841771493, 0.00020249813679, -0.00014303345233, 0.00001406490901,
|
|
|
+ 0.00023699452868, 0.00008661757602, 0.00025744654704, -0.00022715188970, -0.00076146807987, 0.00055185536621, -0.00012325309217, -0.00072356045712, -0.00160693109501, 0.00246682553552, -0.14175094664097, -0.36207047104836, -0.04089594259858, 0.00060774467420, 0.00088646943914, 0.00004865296432, -0.00041878610500, -0.00023025234987, -0.00009676301852, -0.00000000000000, 0.00008409228758,
|
|
|
+ 0.00011432896281, -0.00000707848403, 0.00004698805787, -0.00043642931269, 0.00081384339137, -0.00065635429928, -0.00011831733718, 0.00017413357273, 0.00224463525228, 0.00478497287259, 0.03294761106372, 0.01078986655921, 0.10731782764196, 0.00075034319889, -0.00009241879889, 0.00055023463210, 0.00006596000458, 0.00005045382932, 0.00014874986664, 0.00000000000000, -0.00015369028552,
|
|
|
+ 0.00001037383754, 0.00009250180301, 0.00026204055757, 0.00007424291834, -0.00047751804232, 0.00029184055165, 0.00050921301590, -0.00004825839278, -0.00029933769838, 0.00279659987427, 0.00210463814437, -0.00618590926751, -0.02400829829276, -0.02316811867058, -0.00086368201301, -0.00032258985448, -0.00018304496189, 0.00008438774967, -0.00008305341908, 0.00000000000000, 0.00013047417451,
|
|
|
+ -0.00001376930322, -0.00001723831701, -0.00011543079017, -0.00022646733851, 0.00013467084500, -0.00004661652201, -0.00008419520600, 0.00035772417323, -0.00011815709877, 0.00028718306567, 0.00092207465786, -0.00317224999890, 0.00061770365573, 0.01017294172198, 0.00294739892706, 0.00014669894881, 0.00015702951350, 0.00003432080121, -0.00008555022214, -0.00000000000000, 0.00000454909878,
|
|
|
+ -0.00000196001542, -0.00003198397462, -0.00004425687075, -0.00004129848094, -0.00003789070615, -0.00027583551127, 0.00025874207495, -0.00002334945384, -0.00007259396807, -0.00008295358566, 0.00011360697681, -0.00101968157105, 0.00046784928418, -0.00208410434425, -0.00313158822246, -0.00046005158219, -0.00010552268213, -0.00005850767775, 0.00003971093611, 0.00000000000000, -0.00005275657168,
|
|
|
+ -0.00001065901233, -0.00001934838656, -0.00001220186732, -0.00002060524639, -0.00000225423423, -0.00001894621164, -0.00001533334580, -0.00001791087379, 0.00008156246622, -0.00008441298269, 0.00021060956351, -0.00030303673702, 0.00075949780876, -0.00010539998038, 0.00109045265708, 0.00068949378328, 0.00009268362192, 0.00003471063246, 0.00001204656473, -0.00000000000000, 0.00001500743110,
|
|
|
+ 0.00000105878155, -0.00000910870767, -0.00000172467264, -0.00000722095228, 0.00000699280463, -0.00002061720625, -0.00000889817693, -0.00001993474507, 0.00000370749740, -0.00000090311920, 0.00002677819793, 0.00043428712524, 0.00210293265991, 0.00018200518389, -0.00009621794743, -0.00035250501242, -0.00012996385340, -0.00002185157609, -0.00001116586463, -0.00000000000000, -0.00000451994811,
|
|
|
+ 0.00000424055270, -0.00000463139304, 0.00000301006116, -0.00000123974939, 0.00000632465435, -0.00002090823000, 0.00001773388794, 0.00000121050368, 0.00001886057362, -0.00001043497195, -0.00002269273500, -0.00021979617304, -0.00001043962493, -0.00116343051195, -0.00004193381756, 0.00007944958634, 0.00007301353617, 0.00002082651736, -0.00000119863023, -0.00000000000000, -0.00001440504820,
|
|
|
+ -0.00000391270805, -0.00000490489265, -0.00000504441778, -0.00000904507579, -0.00000111389932, 0.00000597532107, 0.00000047090245, -0.00001553130096, -0.00001524566323, -0.00000522222899, -0.00007707672921, -0.00004165665086, 0.00015764687851, 0.00035649110214, 0.00038701237645, 0.00002386798405, -0.00001946414341, -0.00000913835174, -0.00000489907188, 0.00000000000000, 0.00000172327657,
|
|
|
+ -0.00000015388650, -0.00000603232729, -0.00000397650865, 0.00000280493782, 0.00000463132073, -0.00000788678426, -0.00000471605335, -0.00000283715985, -0.00000422824724, 0.00000366817630, -0.00001159603562, -0.00001625759251, 0.00049116823357, 0.00005048640014, -0.00020234247495, -0.00006341376866, -0.00000807822744, 0.00000070463199, 0.00000014041755, 0.00000000000000, -0.00000718306910};
|
|
|
+#else
|
|
|
+ Real Aspect_ratio = 5;
|
|
|
+ Real coeffmat[21][21] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, Aspect_ratio, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0,
|
|
|
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0};
|
|
|
+#endif
|
|
|
+ Z = 0;
|
|
|
+ Real R = 0;
|
|
|
+ for (long i = -10; i <= 10; i++) {
|
|
|
+ for (long j = -10; j <= 10; j++) {
|
|
|
+ R += coeffmat[i+10][j+10] * sctl::cos(-i*phi + Nperiod*j*theta);
|
|
|
+ Z += coeffmat[i+10][j+10] * sctl::sin(-i*phi + Nperiod*j*theta);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ X = R * sctl::cos(theta);
|
|
|
+ Y = R * sctl::sin(theta);
|
|
|
+ }
|
|
|
+
|
|
|
+ GenericKernel<BiotSavart3D > BiotSavart ;
|
|
|
+ GenericKernel<Laplace3D_FxU > Laplace_FxU ;
|
|
|
+ GenericKernel<Laplace3D_DxU > Laplace_DxU ;
|
|
|
+ GenericKernel<Laplace3D_FxdU> Laplace_FxdU;
|
|
|
+ GenericKernel<Laplace3D_dUxF> Laplace_dUxF;
|
|
|
+ GenericKernel<Laplace3D_Fxd2U> Laplace_Fxd2U;
|
|
|
+ GenericKernel<Laplace3D_dUxD> Laplace_dUxD;
|
|
|
+ GenericKernel<Laplace3D_DxdU> Laplace_DxdU;
|
|
|
+ Quadrature<Real> quadrature_FxU ;
|
|
|
+ Quadrature<Real> quadrature_DxU ;
|
|
|
+ Quadrature<Real> quadrature_FxdU;
|
|
|
+ Quadrature<Real> quadrature_dUxF;
|
|
|
+ Quadrature<Real> quadrature_Fxd2U;
|
|
|
+ Quadrature<Real> quadrature_dUxD;
|
|
|
+
|
|
|
+ ElemLst elements;
|
|
|
+ Vector<Long> NtNp_;
|
|
|
+ Vector<Long> elem_dsp;
|
|
|
+};
|
|
|
+
|
|
|
+template <class Real, Integer ORDER=5> class Spheres {
|
|
|
+ static constexpr Integer COORD_DIM = 3;
|
|
|
+ static constexpr Integer ELEM_DIM = COORD_DIM-1;
|
|
|
+ using PotentialBasis = Basis<Real, ELEM_DIM, ORDER>;
|
|
|
+ using DensityBasis = Basis<Real, ELEM_DIM, ORDER>;
|
|
|
+ using CoordBasis = Basis<Real, ELEM_DIM, ORDER>;
|
|
|
+ using ElemLst = ElemList<COORD_DIM, CoordBasis>;
|
|
|
+
|
|
|
+ public:
|
|
|
+ Spheres(Long N = 0) {
|
|
|
+ Vector<Real> X(N*COORD_DIM);
|
|
|
+ Vector<Real> R(N);
|
|
|
+ X=0;
|
|
|
+ R=1;
|
|
|
+ for (Long i = 0; i < N; i++) X[i*COORD_DIM] = (i==0?-1.015:1.015); ///////////
|
|
|
+ InitSpheres(X,R);
|
|
|
+ }
|
|
|
+
|
|
|
+ const ElemLst& GetElem() const {
|
|
|
+ return elements;
|
|
|
+ }
|
|
|
+
|
|
|
+ static void test() {
|
|
|
+ constexpr Integer order_singular = 35;
|
|
|
+ constexpr Integer order_direct = 35;
|
|
|
+ Comm comm = Comm::World();
|
|
|
+ Profile::Enable(true);
|
|
|
+
|
|
|
+ Long Ns = 2;
|
|
|
+ Spheres S(Ns);
|
|
|
+ S.quadrature_FxT.template Setup<DensityBasis, PotentialBasis>(S.GetElem(), S.Stokes_FxT, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_FxU.template Setup<DensityBasis, PotentialBasis>(S.GetElem(), S.Stokes_FxU, order_singular, order_direct, -1.0, comm);
|
|
|
+ S.quadrature_DxU.template Setup<DensityBasis, PotentialBasis>(S.GetElem(), S.Stokes_DxU, order_singular, order_direct, -1.0, comm);
|
|
|
+
|
|
|
+ const auto SetMotion = [&S](Vector<DensityBasis>& density, const Vector<Real>& force_avg, const Vector<Real>& torque_avg) {
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+ Long Nsurf = S.elem_cnt.Dim();
|
|
|
+ const auto& X = S.GetElem().ElemVector();
|
|
|
+
|
|
|
+ Vector<Real> area, Xc;
|
|
|
+ Vector<DensityBasis> one(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < DensityBasis::Size(); j++) {
|
|
|
+ one[i][j] = 1;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ S.SurfInteg(area, one);
|
|
|
+ S.SurfInteg(Xc, S.GetElem().ElemVector());
|
|
|
+ for (Long i = 0; i < Nsurf; i++) {
|
|
|
+ for (Long k = 0; k < COORD_DIM; k++) {
|
|
|
+ Xc[i*COORD_DIM+k] /= area[i];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ if (density.Dim() != Nelem*COORD_DIM) density.ReInit(Nelem*COORD_DIM);
|
|
|
+ Long elem_itr = 0;
|
|
|
+ for (Long i = 0; i < Nsurf; i++) {
|
|
|
+ for (Long j = 0; j < S.elem_cnt[i]; j++) {
|
|
|
+ for (Long k = 0; k < DensityBasis::Size(); k++) {
|
|
|
+ StaticArray<Real,COORD_DIM> dX;
|
|
|
+ dX[0] = (X[elem_itr*COORD_DIM+0][k] - Xc[i*COORD_DIM+0]);
|
|
|
+ dX[1] = (X[elem_itr*COORD_DIM+1][k] - Xc[i*COORD_DIM+1]);
|
|
|
+ dX[2] = (X[elem_itr*COORD_DIM+2][k] - Xc[i*COORD_DIM+2]);
|
|
|
+ density[elem_itr*COORD_DIM+0][k] = force_avg[i*COORD_DIM+0]*(1/area[i]) + (torque_avg[i*COORD_DIM+1] * dX[2] - torque_avg[i*COORD_DIM+2] * dX[1]) / (2*area[i]/3);
|
|
|
+ density[elem_itr*COORD_DIM+1][k] = force_avg[i*COORD_DIM+1]*(1/area[i]) + (torque_avg[i*COORD_DIM+2] * dX[0] - torque_avg[i*COORD_DIM+0] * dX[2]) / (2*area[i]/3);
|
|
|
+ density[elem_itr*COORD_DIM+2][k] = force_avg[i*COORD_DIM+2]*(1/area[i]) + (torque_avg[i*COORD_DIM+0] * dX[1] - torque_avg[i*COORD_DIM+1] * dX[0]) / (2*area[i]/3);
|
|
|
+ }
|
|
|
+ elem_itr++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ const auto GetMotion = [&S](Vector<Real>& force_avg, Vector<Real>& torque_avg, const Vector<DensityBasis>& density) {
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+ Long Nsurf = S.elem_cnt.Dim();
|
|
|
+ const auto& X = S.GetElem().ElemVector();
|
|
|
+
|
|
|
+ S.SurfInteg(force_avg, density);
|
|
|
+
|
|
|
+ Vector<Real> area, Xc;
|
|
|
+ Vector<DensityBasis> one(Nelem);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < DensityBasis::Size(); j++) {
|
|
|
+ one[i][j] = 1;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ S.SurfInteg(area, one);
|
|
|
+ S.SurfInteg(Xc, S.GetElem().ElemVector());
|
|
|
+ for (Long i = 0; i < Nsurf; i++) {
|
|
|
+ for (Long k = 0; k < COORD_DIM; k++) {
|
|
|
+ Xc[i*COORD_DIM+k] /= area[i];
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ { // Set torque_avg
|
|
|
+ Long elem_itr = 0;
|
|
|
+ Vector<DensityBasis> torque(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nsurf; i++) {
|
|
|
+ for (Long j = 0; j < S.elem_cnt[i]; j++) {
|
|
|
+ for (Long k = 0; k < DensityBasis::Size(); k++) {
|
|
|
+ StaticArray<Real,COORD_DIM> dX;
|
|
|
+ dX[0] = (X[elem_itr*COORD_DIM+0][k] - Xc[i*COORD_DIM+0]);
|
|
|
+ dX[1] = (X[elem_itr*COORD_DIM+1][k] - Xc[i*COORD_DIM+1]);
|
|
|
+ dX[2] = (X[elem_itr*COORD_DIM+2][k] - Xc[i*COORD_DIM+2]);
|
|
|
+ torque[elem_itr*COORD_DIM+0][k] = dX[1] * density[elem_itr*COORD_DIM+2][k] - dX[2] * density[elem_itr*COORD_DIM+1][k];
|
|
|
+ torque[elem_itr*COORD_DIM+1][k] = dX[2] * density[elem_itr*COORD_DIM+0][k] - dX[0] * density[elem_itr*COORD_DIM+2][k];
|
|
|
+ torque[elem_itr*COORD_DIM+2][k] = dX[0] * density[elem_itr*COORD_DIM+1][k] - dX[1] * density[elem_itr*COORD_DIM+0][k];
|
|
|
+ }
|
|
|
+ elem_itr++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ S.SurfInteg(torque_avg, torque);
|
|
|
+ }
|
|
|
+ };
|
|
|
+ const auto BIOpL = [&GetMotion,&SetMotion](Vector<DensityBasis>& potential, const Vector<DensityBasis>& density) {
|
|
|
+ Vector<Real> force_avg, torque_avg;
|
|
|
+ GetMotion(force_avg, torque_avg, density);
|
|
|
+ SetMotion(potential, force_avg, torque_avg);
|
|
|
+ };
|
|
|
+ const auto BIOpK = [&S](Vector<DensityBasis>& potential, const Vector<DensityBasis>& density) {
|
|
|
+ Vector<DensityBasis> traction;
|
|
|
+ S.quadrature_FxT.Eval(traction, S.GetElem(), density, S.Stokes_FxT);
|
|
|
+
|
|
|
+ Vector<CoordBasis> dX;
|
|
|
+ const auto X = S.GetElem().ElemVector();
|
|
|
+ CoordBasis::Grad(dX, X);
|
|
|
+
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+ Long Nnodes = CoordBasis::Size();
|
|
|
+ potential.ReInit(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ StaticArray<Real,COORD_DIM> Xn;
|
|
|
+ Xn[0] = dX[i*COORD_DIM*2+2][j]*dX[i*COORD_DIM*2+5][j] - dX[i*COORD_DIM*2+4][j]*dX[i*COORD_DIM*2+3][j];
|
|
|
+ Xn[1] = dX[i*COORD_DIM*2+4][j]*dX[i*COORD_DIM*2+1][j] - dX[i*COORD_DIM*2+0][j]*dX[i*COORD_DIM*2+5][j];
|
|
|
+ Xn[2] = dX[i*COORD_DIM*2+0][j]*dX[i*COORD_DIM*2+3][j] - dX[i*COORD_DIM*2+2][j]*dX[i*COORD_DIM*2+1][j];
|
|
|
+ Real AreaElem = sqrt<Real>(Xn[0]*Xn[0] + Xn[1]*Xn[1] + Xn[2]*Xn[2]);
|
|
|
+ Real OOAreaElem = 1 / AreaElem;
|
|
|
+ Xn[0] *= OOAreaElem;
|
|
|
+ Xn[1] *= OOAreaElem;
|
|
|
+ Xn[2] *= OOAreaElem;
|
|
|
+
|
|
|
+ potential[i*COORD_DIM+0][j] = traction[i*COORD_DIM*COORD_DIM+0][j]*Xn[0] + traction[i*COORD_DIM*COORD_DIM+1][j]*Xn[1] + traction[i*COORD_DIM*COORD_DIM+2][j]*Xn[2];
|
|
|
+ potential[i*COORD_DIM+1][j] = traction[i*COORD_DIM*COORD_DIM+3][j]*Xn[0] + traction[i*COORD_DIM*COORD_DIM+4][j]*Xn[1] + traction[i*COORD_DIM*COORD_DIM+5][j]*Xn[2];
|
|
|
+ potential[i*COORD_DIM+2][j] = traction[i*COORD_DIM*COORD_DIM+6][j]*Xn[0] + traction[i*COORD_DIM*COORD_DIM+7][j]*Xn[1] + traction[i*COORD_DIM*COORD_DIM+8][j]*Xn[2];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ const auto BIOp_half_K_L = [&S,&BIOpK,&BIOpL](Vector<DensityBasis>& potential, const Vector<DensityBasis>& density) {
|
|
|
+ Vector<DensityBasis> potential_K;
|
|
|
+ Vector<DensityBasis> potential_L;
|
|
|
+ BIOpK(potential_K, density);
|
|
|
+ BIOpL(potential_L, density);
|
|
|
+
|
|
|
+ if (potential.Dim() != potential_K.Dim()) {
|
|
|
+ potential.ReInit(potential_K.Dim());
|
|
|
+ }
|
|
|
+ for (Long i = 0; i < potential_K.Dim(); i++) {
|
|
|
+ for (Long k = 0; k < DensityBasis::Size(); k++) {
|
|
|
+ potential[i][k] = -0.5*density[i][k] + potential_K[i][k] + potential_L[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ const auto BIOp_half_K = [&S,&BIOpK,&BIOpL](Vector<DensityBasis>& potential, const Vector<DensityBasis>& density) {
|
|
|
+ Vector<DensityBasis> potential_K;
|
|
|
+ BIOpK(potential_K, density);
|
|
|
+
|
|
|
+ if (potential.Dim() != potential_K.Dim()) {
|
|
|
+ potential.ReInit(potential_K.Dim());
|
|
|
+ }
|
|
|
+ for (Long i = 0; i < potential_K.Dim(); i++) {
|
|
|
+ for (Long k = 0; k < DensityBasis::Size(); k++) {
|
|
|
+ potential[i][k] = -0.5*density[i][k] + potential_K[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+ const auto BIOp_half_S_D = [&S,&BIOpL](Vector<DensityBasis>& potential, const Vector<DensityBasis>& density) {
|
|
|
+ Vector<DensityBasis> U;
|
|
|
+ S.quadrature_DxU.Eval(U, S.GetElem(), density, S.Stokes_DxU);
|
|
|
+
|
|
|
+ Vector<PotentialBasis> U1;
|
|
|
+ Vector<DensityBasis> sigma1;
|
|
|
+ BIOpL(sigma1,density);
|
|
|
+ S.quadrature_FxU.Eval(U1, S.GetElem(), sigma1, S.Stokes_FxU);
|
|
|
+
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+ Long Nnodes = CoordBasis::Size();
|
|
|
+ potential.ReInit(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ potential[i*COORD_DIM+0][j] = 0.5*density[i*COORD_DIM+0][j] + U[i*COORD_DIM+0][j] + U1[i*COORD_DIM+0][j];
|
|
|
+ potential[i*COORD_DIM+1][j] = 0.5*density[i*COORD_DIM+1][j] + U[i*COORD_DIM+1][j] + U1[i*COORD_DIM+1][j];
|
|
|
+ potential[i*COORD_DIM+2][j] = 0.5*density[i*COORD_DIM+2][j] + U[i*COORD_DIM+2][j] + U1[i*COORD_DIM+2][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ Vector<PotentialBasis> U;
|
|
|
+ { // Rachh
|
|
|
+ Vector<DensityBasis> sigma0;
|
|
|
+ { // Set sigma0
|
|
|
+ srand48(comm.Rank());
|
|
|
+ Vector<Real> force(Ns*COORD_DIM), torque(Ns*COORD_DIM);
|
|
|
+ //for (auto& x : force) x = drand48();
|
|
|
+ //for (auto& x : torque) x = drand48();
|
|
|
+ force = 0;
|
|
|
+ torque = 0;
|
|
|
+ force[0] = 1;
|
|
|
+ //force[4] = 1;
|
|
|
+ SetMotion(sigma0, force, torque);
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<DensityBasis> rhs;
|
|
|
+ BIOp_half_K(rhs, sigma0);
|
|
|
+
|
|
|
+ Vector<DensityBasis> sigma;
|
|
|
+ { // Set sigma
|
|
|
+ Long Nnode = DensityBasis::Size();
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp A = [&S,&BIOp_half_K_L](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
|
|
|
+ Long Nnode = DensityBasis::Size();
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+ Ax->ReInit(Nelem*COORD_DIM*Nnode);
|
|
|
+
|
|
|
+ Vector<DensityBasis> x_(Nelem*COORD_DIM), Ax_(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) { // Set x_
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ x_[i][k] = x[i*Nnode+k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ BIOp_half_K_L(Ax_, x_);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) { // Set Ax
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ (*Ax)[i*Nnode+k] = Ax_[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ Vector<Real> sigma_(Nelem*COORD_DIM*Nnode), rhs_(Nelem*COORD_DIM*Nnode);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) {// Set rhs_
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ rhs_[i*Nnode+k] = rhs[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ sigma_ = 0;
|
|
|
+
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&sigma_, A, rhs_, 1e-6, 50);
|
|
|
+ sigma.ReInit(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) {// Set sigma
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ sigma[i][k] = sigma_[i*Nnode+k] - sigma0[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ S.quadrature_FxU.Eval(U, S.GetElem(), sigma, S.Stokes_FxU);
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu_sigma;
|
|
|
+ vtu_sigma.AddElems(S.elements, sigma, ORDER);
|
|
|
+ vtu_sigma.WriteVTK("sphere-sigma0", comm);
|
|
|
+
|
|
|
+ VTUData vtu_U;
|
|
|
+ vtu_U.AddElems(S.elements, U, ORDER);
|
|
|
+ vtu_U.WriteVTK("sphere-U0", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ { // Tornberg
|
|
|
+ Vector<DensityBasis> rhs;
|
|
|
+ BIOpL(rhs, U);
|
|
|
+
|
|
|
+ Vector<DensityBasis> sigma;
|
|
|
+ { // Set sigma
|
|
|
+ Long Nnode = DensityBasis::Size();
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+
|
|
|
+ typename sctl::ParallelSolver<Real>::ParallelOp A = [&S,&BIOp_half_S_D](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
|
|
|
+ Long Nnode = DensityBasis::Size();
|
|
|
+ Long Nelem = S.GetElem().NElem();
|
|
|
+ Ax->ReInit(Nelem*COORD_DIM*Nnode);
|
|
|
+
|
|
|
+ Vector<DensityBasis> x_(Nelem*COORD_DIM), Ax_(Nelem*COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) { // Set x_
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ x_[i][k] = x[i*Nnode+k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ BIOp_half_S_D(Ax_, x_);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) { // Set Ax
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ (*Ax)[i*Nnode+k] = Ax_[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ };
|
|
|
+
|
|
|
+ Vector<Real> sigma_(Nelem*COORD_DIM*Nnode), rhs_(Nelem*COORD_DIM*Nnode);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) {// Set rhs_
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ rhs_[i*Nnode+k] = rhs[i][k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ sigma_ = 0;
|
|
|
+
|
|
|
+ ParallelSolver<Real> linear_solver(comm, true);
|
|
|
+ linear_solver(&sigma_, A, rhs_, 1e-6, 50);
|
|
|
+ sigma.ReInit(Nelem * COORD_DIM);
|
|
|
+ for (Long i = 0; i < Nelem*COORD_DIM; i++) {// Set sigma
|
|
|
+ for (Long k = 0; k < Nnode; k++) {
|
|
|
+ sigma[i][k] = sigma_[i*Nnode+k];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Vector<PotentialBasis> U1;
|
|
|
+ BIOp_half_S_D(U1, sigma);
|
|
|
+ { // Write VTU
|
|
|
+ VTUData vtu_sigma;
|
|
|
+ vtu_sigma.AddElems(S.elements, sigma, ORDER);
|
|
|
+ vtu_sigma.WriteVTK("sphere-sigma1", comm);
|
|
|
+
|
|
|
+ VTUData vtu_U;
|
|
|
+ vtu_U.AddElems(S.elements, U1, ORDER);
|
|
|
+ vtu_U.WriteVTK("sphere-U1", comm);
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Profile::print(&comm);
|
|
|
+ }
|
|
|
+
|
|
|
+ private:
|
|
|
+
|
|
|
+ template <class FnBasis> void SurfInteg(Vector<Real>& I, const Vector<FnBasis>& f) {
|
|
|
+ static_assert(std::is_same<FnBasis,CoordBasis>::value, "FnBasis is different from CoordBasis");
|
|
|
+ const Long Nelem = elements.NElem();
|
|
|
+ const Long dof = f.Dim() / Nelem;
|
|
|
+ SCTL_ASSERT(f.Dim() == Nelem * dof);
|
|
|
+
|
|
|
+ auto nodes = FnBasis::Nodes();
|
|
|
+ auto quad_wts = FnBasis::QuadWts();
|
|
|
+ const Long Nnodes = FnBasis::Size();
|
|
|
+ auto EvalOp = CoordBasis::SetupEval(nodes);
|
|
|
+
|
|
|
+ Vector<CoordBasis> dX;
|
|
|
+ const auto& X = elements.ElemVector();
|
|
|
+ SCTL_ASSERT(X.Dim() == Nelem * COORD_DIM);
|
|
|
+ CoordBasis::Grad(dX, X);
|
|
|
+
|
|
|
+ Matrix<Real> I_(Nelem, dof);
|
|
|
+ for (Long i = 0; i < Nelem; i++) {
|
|
|
+ for (Long k = 0; k < dof; k++) {
|
|
|
+ I_[i][k] = 0;
|
|
|
+ }
|
|
|
+ for (Long j = 0; j < Nnodes; j++) {
|
|
|
+ Real dA = 0;
|
|
|
+ StaticArray<Real,COORD_DIM> Xn;
|
|
|
+ Xn[0] = dX[i*COORD_DIM*2+2][j] * dX[i*COORD_DIM*2+5][j] - dX[i*COORD_DIM*2+3][j] * dX[i*COORD_DIM*2+4][j];
|
|
|
+ Xn[1] = dX[i*COORD_DIM*2+4][j] * dX[i*COORD_DIM*2+1][j] - dX[i*COORD_DIM*2+5][j] * dX[i*COORD_DIM*2+0][j];
|
|
|
+ Xn[2] = dX[i*COORD_DIM*2+0][j] * dX[i*COORD_DIM*2+3][j] - dX[i*COORD_DIM*2+1][j] * dX[i*COORD_DIM*2+2][j];
|
|
|
+ dA += sqrt<Real>(Xn[0]*Xn[0] + Xn[1]*Xn[1] + Xn[2]*Xn[2]) * quad_wts[j];
|
|
|
+ for (Long k = 0; k < dof; k++) {
|
|
|
+ I_[i][k] += dA * f[i*dof+k][j];
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ Long Ns = elem_cnt.Dim();
|
|
|
+ if (I.Dim() != Ns * dof) I.ReInit(Ns * dof);
|
|
|
+ I = 0;
|
|
|
+ Long elem_itr = 0;
|
|
|
+ for (Long i = 0; i < Ns; i++) {
|
|
|
+ for (Long j = 0; j < elem_cnt[i]; j++) {
|
|
|
+ for (Long k = 0; k < dof; k++) {
|
|
|
+ I[i*dof+k] += I_[elem_itr][k];
|
|
|
+ }
|
|
|
+ elem_itr++;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ void InitSpheres(const Vector<Real> X, const Vector<Real>& R){
|
|
|
+ SCTL_ASSERT(X.Dim() == R.Dim() * COORD_DIM);
|
|
|
+ Long N = R.Dim();
|
|
|
+ elements.ReInit(2*COORD_DIM*N);
|
|
|
+ auto nodes = ElemLst::CoordBasis::Nodes();
|
|
|
+ for (Long l = 0; l < N; l++) {
|
|
|
+ for (Integer i = 0; i < COORD_DIM; i++) {
|
|
|
+ for (Integer j = 0; j < 2; j++) {
|
|
|
+ for (int k = 0; k < ElemLst::CoordBasis::Size(); k++) {
|
|
|
+ Real coord[COORD_DIM];
|
|
|
+ coord[(i+0)%COORD_DIM] = (j ? -1.0 : 1.0);
|
|
|
+ coord[(i+1)%COORD_DIM] = 2.0 * nodes[j?1:0][k] - 1.0;
|
|
|
+ coord[(i+2)%COORD_DIM] = 2.0 * nodes[j?0:1][k] - 1.0;
|
|
|
+ Real R0 = sqrt<Real>(coord[0]*coord[0] + coord[1]*coord[1] + coord[2]*coord[2]);
|
|
|
+
|
|
|
+ elements((l*COORD_DIM+i)*2+j,0)[k] = X[l*COORD_DIM+0] + R[l] * coord[0] / R0;
|
|
|
+ elements((l*COORD_DIM+i)*2+j,1)[k] = X[l*COORD_DIM+1] + R[l] * coord[1] / R0;
|
|
|
+ elements((l*COORD_DIM+i)*2+j,2)[k] = X[l*COORD_DIM+2] + R[l] * coord[2] / R0;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ elem_cnt.ReInit(N);
|
|
|
+ elem_cnt = 6;
|
|
|
+ }
|
|
|
+
|
|
|
+ GenericKernel<Stokes3D_DxU> Stokes_DxU;
|
|
|
+ GenericKernel<Stokes3D_FxU> Stokes_FxU;
|
|
|
+ GenericKernel<Stokes3D_FxT> Stokes_FxT;
|
|
|
+
|
|
|
+ Quadrature<Real> quadrature_DxU;
|
|
|
+ Quadrature<Real> quadrature_FxU;
|
|
|
+ Quadrature<Real> quadrature_FxT;
|
|
|
+
|
|
|
+ ElemLst elements;
|
|
|
+ Vector<Long> elem_cnt;
|
|
|
+};
|
|
|
+
|
|
|
} // end namespace
|
|
|
|
|
|
#endif //_SCTL_BOUNDARY_QUADRATURE_HPP_
|
