PolyhedralFace Class Reference

Collaboration diagram for PolyhedralFace:

Public Member Functions
	PolyhedralFace (const std::vector< kvector_t > &_V=std::vector< kvector_t >(), bool _sym_S2=false)
double	area () const
double	pyramidalVolume () const
double	radius3d () const
complex_t	normalProjectionConj (cvector_t q) const
complex_t	ff_n (int m, cvector_t q) const
complex_t	ff (cvector_t q, bool sym_Ci) const
complex_t	ff_2D (cvector_t qpa) const
void	assert_Ci (const PolyhedralFace &other) const

Static Public Member Functions
static double	diameter (const std::vector< kvector_t > &V)

Private Member Functions
void	decompose_q (cvector_t q, complex_t &qperp, cvector_t &qpa) const
complex_t	ff_n_core (int m, cvector_t qpa, complex_t qperp) const
complex_t	edge_sum_ff (cvector_t q, cvector_t qpa, bool sym_Ci) const
complex_t	expansion (complex_t fac_even, complex_t fac_odd, cvector_t qpa, double abslevel) const

Private Attributes
bool	sym_S2
std::vector< PolyhedralEdge >	edges
double	m_area
kvector_t	m_normal
double	m_rperp
double	m_radius_2d
double	m_radius_3d

Static Private Attributes
static double	qpa_limit_series = 3e-2
static int	n_limit_series = 20

Detailed Description

A polygon, for form factor computation.

Definition at line 43 of file PolyhedralComponents.h.

Constructor & Destructor Documentation

PolyhedralFace()

PolyhedralFace::PolyhedralFace	(	const std::vector< kvector_t > &	V = std::vector<kvector_t>(),
		bool	_sym_S2 = false
	)

Sets internal variables for given vertex chain.

Parameters

V	oriented vertex list
_sym_S2	true if face has a perpedicular two-fold symmetry axis

Definition at line 120 of file PolyhedralComponents.cpp.

                                                                          : sym_S2(_sym_S2)
{
    size_t NV = V.size();
    if (!NV)
        throw std::logic_error("Face with no edges");
    if (NV < 3)
        throw std::logic_error("Face with less than three edges");
 
    // compute radius in 2d and 3d
    m_radius_2d = diameter(V) / 2;
    m_radius_3d = 0;
    for (const kvector_t& v : V)
        m_radius_3d = std::max(m_radius_3d, v.mag());
 
    // Initialize list of 'edges'.
    // Do not create an edge if two vertices are too close to each other.
    // TODO This is implemented in a somewhat sloppy way: we just skip an edge if it would
    //      be too short. This leaves tiny open edges. In a clean implementation, we
    //      rather should merge adjacent vertices before generating edges.
    for (size_t j = 0; j < NV; ++j) {
        size_t jj = (j + 1) % NV;
        if ((V[j] - V[jj]).mag() < 1e-14 * m_radius_2d)
            continue; // distance too short -> skip this edge
        edges.push_back(PolyhedralEdge(V[j], V[jj]));
    }
    size_t NE = edges.size();
    if (NE < 3)
        throw std::invalid_argument("Face has less than three non-vanishing edges");
 
    // compute n_k, rperp
    m_normal = kvector_t();
    for (size_t j = 0; j < NE; ++j) {
        size_t jj = (j + 1) % NE;
        kvector_t ee = edges[j].E().cross(edges[jj].E());
        if (ee.mag2() == 0) {
            throw std::logic_error("Two adjacent edges are parallel");
        }
        m_normal += ee.unit();
    }
    m_normal /= NE;
    m_rperp = 0;
    for (size_t j = 0; j < NV; ++j)
        m_rperp += V[j].dot(m_normal);
    m_rperp /= NV;
    // assert that the vertices lay in a plane
    for (size_t j = 1; j < NV; ++j)
        if (std::abs(V[j].dot(m_normal) - m_rperp) > 1e-14 * m_radius_3d)
            throw std::logic_error("Face is not planar");
    // compute m_area
    m_area = 0;
    for (size_t j = 0; j < NV; ++j) {
        size_t jj = (j + 1) % NV;
        m_area += m_normal.dot(V[j].cross(V[jj])) / 2;
    }
    // only now deal with inversion symmetry
    if (sym_S2) {
        if (NE & 1)
            throw std::logic_error("Odd #edges violates symmetry S2");
        NE /= 2;
        for (size_t j = 0; j < NE; ++j) {
            if (((edges[j].R() - m_rperp * m_normal) + (edges[j + NE].R() - m_rperp * m_normal))
                    .mag()
                > 1e-12 * m_radius_2d)
                throw std::logic_error("Edge centers violate symmetry S2");
            if ((edges[j].E() + edges[j + NE].E()).mag() > 1e-12 * m_radius_2d)
                throw std::logic_error("Edge vectors violate symmetry S2");
        }
        // keep only half of the egdes
        edges.erase(edges.begin() + NE, edges.end());
    }
}

References diameter(), BasicVector3D< T >::dot(), edges, m_area, m_normal, m_radius_2d, m_radius_3d, m_rperp, BasicVector3D< T >::mag2(), sym_S2, and BasicVector3D< T >::unit().

Here is the call graph for this function:

Member Function Documentation

diameter()

double PolyhedralFace::diameter ( const std::vector< kvector_t > &  V )

static

Static method, returns diameter of circle that contains all vertices.

Definition at line 98 of file PolyhedralComponents.cpp.

{
    double diameterFace = 0;
    for (size_t j = 0; j < V.size(); ++j)
        for (size_t jj = j + 1; jj < V.size(); ++jj)
            diameterFace = std::max(diameterFace, (V[j] - V[jj]).mag());
    return diameterFace;
}

Referenced by PolyhedralFace(), and Polyhedron::Polyhedron().

area()

double PolyhedralFace::area ( ) const

inline

Definition at line 54 of file PolyhedralComponents.h.

{ return m_area; }

References m_area.

pyramidalVolume()

double PolyhedralFace::pyramidalVolume ( ) const

inline

Definition at line 55 of file PolyhedralComponents.h.

{ return m_rperp * m_area / 3; }

References m_area, and m_rperp.

radius3d()

double PolyhedralFace::radius3d ( ) const

inline

Definition at line 56 of file PolyhedralComponents.h.

{ return m_radius_3d; }

References m_radius_3d.

normalProjectionConj()

complex_t PolyhedralFace::normalProjectionConj ( cvector_t  q ) const

inline

Returns conj(q)*normal [BasicVector3D::dot is antilinear in 'this' argument].

Definition at line 58 of file PolyhedralComponents.h.

{ return q.dot(m_normal); }

References BasicVector3D< T >::dot(), and m_normal.

Here is the call graph for this function:

ff_n()

complex_t PolyhedralFace::ff_n	(	int	m,
		cvector_t	q
	)		const

Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor.

Definition at line 235 of file PolyhedralComponents.cpp.

{
    complex_t qn = q.dot(m_normal); // conj(q)*normal (dot is antilinear in 'this' argument)
    if (std::abs(qn) < eps * q.mag())
        return 0.;
    complex_t qperp;
    cvector_t qpa;
    decompose_q(q, qperp, qpa);
    double qpa_mag2 = qpa.mag2();
    if (qpa_mag2 == 0.) {
        return qn * pow(qperp * m_rperp, n) * m_area * ReciprocalFactorialArray[n];
    } else if (sym_S2) {
        return qn * (ff_n_core(n, qpa, qperp) + ff_n_core(n, -qpa, qperp)) / qpa_mag2;
    } else {
        complex_t tmp = ff_n_core(n, qpa, qperp);
#ifdef POLYHEDRAL_DIAGNOSTIC
        if (diagnosis.debmsg >= 3)
            std::cout << "DBX ff_n " << n << " " << qn << " " << tmp << " " << qpa_mag2 << "\n";
#endif
        return qn * tmp / qpa_mag2;
    }
}

References decompose_q(), BasicVector3D< T >::dot(), anonymous_namespace{PolyhedralComponents.cpp}::eps, ff_n_core(), m_area, m_normal, m_rperp, BasicVector3D< T >::mag(), BasicVector3D< T >::mag2(), anonymous_namespace{PolyhedralComponents.cpp}::ReciprocalFactorialArray, and sym_S2.

Here is the call graph for this function:

ff()

complex_t PolyhedralFace::ff	(	cvector_t	q,
		bool	sym_Ci
	)		const

Returns the contribution ff(q) of this face to the polyhedral form factor.

Definition at line 328 of file PolyhedralComponents.cpp.

{
    complex_t qperp;
    cvector_t qpa;
    decompose_q(q, qperp, qpa);
    double qpa_red = m_radius_2d * qpa.mag();
    complex_t qr_perp = qperp * m_rperp;
    complex_t ff0 = (sym_Ci ? 2. * I * sin(qr_perp) : exp_I(qr_perp)) * m_area;
    if (qpa_red == 0) {
        return ff0;
    } else if (qpa_red < qpa_limit_series && !sym_S2) {
        // summation of power series
        complex_t fac_even;
        complex_t fac_odd;
        if (sym_Ci) {
            fac_even = 2. * mul_I(sin(qr_perp));
            fac_odd = 2. * cos(qr_perp);
        } else {
            fac_even = exp_I(qr_perp);
            fac_odd = fac_even;
        }
        return ff0 + expansion(fac_even, fac_odd, qpa, std::abs(ff0));
    } else {
        // direct evaluation of analytic formula
        complex_t prefac;
        if (sym_S2)
            prefac = sym_Ci ? -8. * sin(qr_perp) : 4. * mul_I(exp_I(qr_perp));
        else
            prefac = sym_Ci ? 4. : 2. * exp_I(qr_perp);
#ifdef POLYHEDRAL_DIAGNOSTIC
        if (diagnosis.debmsg >= 2)
            std::cout << "       qrperp=" << qr_perp << " => prefac=" << prefac << "\n";
#endif
        return prefac * edge_sum_ff(q, qpa, sym_Ci) / mul_I(qpa.mag2());
    }
}

References decompose_q(), edge_sum_ff(), exp_I(), expansion(), I, m_area, m_radius_2d, m_rperp, BasicVector3D< T >::mag(), BasicVector3D< T >::mag2(), mul_I(), qpa_limit_series, and sym_S2.

Referenced by ff_2D().

Here is the call graph for this function:

ff_2D()

complex_t PolyhedralFace::ff_2D

(

cvector_t

qpa

)

const

Returns the two-dimensional form factor of this face, for use in a prism.

Definition at line 367 of file PolyhedralComponents.cpp.

{
    if (std::abs(qpa.dot(m_normal)) > eps * qpa.mag())
        throw std::logic_error("ff_2D called with perpendicular q component");
    double qpa_red = m_radius_2d * qpa.mag();
    if (qpa_red == 0) {
        return m_area;
    } else if (qpa_red < qpa_limit_series && !sym_S2) {
        // summation of power series
        return m_area + expansion(1., 1., qpa, std::abs(m_area));
    } else {
        // direct evaluation of analytic formula
        complex_t ff = edge_sum_ff(qpa, qpa, false);
        complex_t ret = (sym_S2 ? 4. : 2. / I) * ff / qpa.mag2();
#ifdef POLYHEDRAL_DIAGNOSTIC
        if (diagnosis.debmsg >= 2)
            std::cout << std::setprecision(16) << "    ret=" << ret << " ff=" << ff << "\n";
#endif
        return ret;
    }
}

References BasicVector3D< T >::dot(), edge_sum_ff(), anonymous_namespace{PolyhedralComponents.cpp}::eps, expansion(), ff(), I, m_area, m_normal, m_radius_2d, BasicVector3D< T >::mag(), BasicVector3D< T >::mag2(), qpa_limit_series, and sym_S2.

Here is the call graph for this function:

assert_Ci()

void PolyhedralFace::assert_Ci

(

const PolyhedralFace &

other

)

const

Throws if deviation from inversion symmetry is detected. Does not check vertices.

Definition at line 391 of file PolyhedralComponents.cpp.

{
    if (std::abs(m_rperp - other.m_rperp) > 1e-15 * (m_rperp + other.m_rperp))
        throw std::logic_error("Faces with different distance from origin violate symmetry Ci");
    if (std::abs(m_area - other.m_area) > 1e-15 * (m_area + other.m_area))
        throw std::logic_error("Faces with different areas violate symmetry Ci");
    if ((m_normal + other.m_normal).mag() > 1e-14)
        throw std::logic_error("Faces do not have opposite orientation, violating symmetry Ci");
}

References m_area, m_normal, and m_rperp.

decompose_q()

void PolyhedralFace::decompose_q ( cvector_t  q, complex_t &  qperp, cvector_t &  qpa  ) const

private

Sets qperp and qpa according to argument q and to this polygon's normal.

Definition at line 194 of file PolyhedralComponents.cpp.

{
    qperp = m_normal.dot(q);
    qpa = q - qperp * m_normal;
    // improve numeric accuracy:
    qpa -= m_normal.dot(qpa) * m_normal;
    if (qpa.mag() < eps * std::abs(qperp))
        qpa = cvector_t(0., 0., 0.);
}

References BasicVector3D< T >::dot(), anonymous_namespace{PolyhedralComponents.cpp}::eps, m_normal, and BasicVector3D< T >::mag().

Referenced by ff(), and ff_n().

Here is the call graph for this function:

ff_n_core()

complex_t PolyhedralFace::ff_n_core ( int  m, cvector_t  qpa, complex_t  qperp  ) const

private

Returns core contribution to f_n.

Definition at line 206 of file PolyhedralComponents.cpp.

{
    cvector_t prevec = 2. * m_normal.cross(qpa); // complex conjugation will take place in .dot
    complex_t ret = 0;
    complex_t vfacsum = 0;
    complex_t qrperp = qperp * m_rperp;
    for (size_t i = 0; i < edges.size(); ++i) {
        const PolyhedralEdge& e = edges[i];
        complex_t vfac;
        if (sym_S2 || i < edges.size() - 1) {
            vfac = prevec.dot(e.E());
            vfacsum += vfac;
        } else {
            vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
        }
        complex_t tmp = e.contrib(m + 1, qpa, qrperp);
        ret += vfac * tmp;
#ifdef POLYHEDRAL_DIAGNOSTIC
        if (diagnosis.debmsg >= 4)
            std::cout << std::scientific << std::showpos << std::setprecision(16)
                      << "DBX ff_n_core " << m << " " << vfac << " " << tmp
                      << " term=" << vfac * tmp << " sum=" << ret << "\n";
#endif
    }
    return ret;
}

References PolyhedralEdge::contrib(), BasicVector3D< T >::cross(), BasicVector3D< T >::dot(), PolyhedralEdge::E(), edges, m_normal, m_rperp, StringUtils::scientific(), and sym_S2.

Referenced by expansion(), and ff_n().

Here is the call graph for this function:

edge_sum_ff()

complex_t PolyhedralFace::edge_sum_ff ( cvector_t  q, cvector_t  qpa, bool  sym_Ci  ) const

private

Returns core contribution to analytic 2d form factor.

Definition at line 296 of file PolyhedralComponents.cpp.

{
    cvector_t prevec = m_normal.cross(qpa); // complex conjugation will take place in .dot
    complex_t sum = 0;
    complex_t vfacsum = 0;
    for (size_t i = 0; i < edges.size(); ++i) {
        const PolyhedralEdge& e = edges[i];
        complex_t qE = e.qE(qpa);
        complex_t qR = e.qR(qpa);
        complex_t Rfac = sym_S2 ? sin(qR) : (sym_Ci ? cos(e.qR(q)) : exp_I(qR));
        complex_t vfac;
        if (sym_S2 || i < edges.size() - 1) {
            vfac = prevec.dot(e.E());
            vfacsum += vfac;
        } else {
            vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
        }
        complex_t term = vfac * MathFunctions::sinc(qE) * Rfac;
        sum += term;
#ifdef POLYHEDRAL_DIAGNOSTIC
        if (diagnosis.debmsg >= 2)
            std::cout << std::scientific << std::showpos << std::setprecision(16)
                      << "    sum=" << sum << " term=" << term << " vf=" << vfac << " qE=" << qE
                      << " qR=" << qR << " sinc=" << MathFunctions::sinc(qE) << " Rfac=" << Rfac
                      << "\n";
#endif
    }
    return sum;
}

References BasicVector3D< T >::cross(), BasicVector3D< T >::dot(), PolyhedralEdge::E(), edges, exp_I(), m_normal, PolyhedralEdge::qE(), PolyhedralEdge::qR(), StringUtils::scientific(), MathFunctions::sinc(), and sym_S2.

Referenced by ff(), and ff_2D().

Here is the call graph for this function:

expansion()

complex_t PolyhedralFace::expansion ( complex_t  fac_even, complex_t  fac_odd, cvector_t  qpa, double  abslevel  ) const

private

Returns sum of n>=1 terms of qpa expansion of 2d form factor.

Definition at line 260 of file PolyhedralComponents.cpp.

{
#ifdef POLYHEDRAL_DIAGNOSTIC
    diagnosis.nExpandedFaces += 1;
#endif
    complex_t sum = 0;
    complex_t n_fac = I;
    int count_return_condition = 0;
    for (int n = 1; n < n_limit_series; ++n) {
#ifdef POLYHEDRAL_DIAGNOSTIC
        diagnosis.maxOrder = std::max(diagnosis.maxOrder, n);
#endif
        complex_t term = n_fac * (n & 1 ? fac_odd : fac_even) * ff_n_core(n, qpa, 0) / qpa.mag2();
#ifdef POLYHEDRAL_DIAGNOSTIC
        if (diagnosis.debmsg >= 2)
            std::cout << std::setprecision(16) << "    sum=" << sum << " +term=" << term << "\n";
#endif
        sum += term;
        if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * abslevel)
            ++count_return_condition;
        else
            count_return_condition = 0;
        if (count_return_condition > 2)
            return sum; // regular exit
        n_fac = mul_I(n_fac);
    }
#ifdef POLYHEDRAL_DIAGNOSTIC
    if (!diagnosis.request_convergence)
        return sum;
#endif
    throw std::runtime_error("Series f(q_pa) not converged");
}

References anonymous_namespace{PolyhedralComponents.cpp}::eps, ff_n_core(), I, BasicVector3D< T >::mag2(), mul_I(), and n_limit_series.

Referenced by ff(), and ff_2D().

Here is the call graph for this function:

Member Data Documentation

qpa_limit_series

double PolyhedralFace::qpa_limit_series = 3e-2

staticprivate

determines when use power series

Definition at line 65 of file PolyhedralComponents.h.

Referenced by ff(), and ff_2D().

n_limit_series

int PolyhedralFace::n_limit_series = 20

staticprivate

Definition at line 66 of file PolyhedralComponents.h.

Referenced by expansion().

sym_S2

bool PolyhedralFace::sym_S2

private

if true, then edges obtainable by inversion are not provided

Definition at line 68 of file PolyhedralComponents.h.

Referenced by edge_sum_ff(), ff(), ff_2D(), ff_n(), ff_n_core(), and PolyhedralFace().

edges

std::vector<PolyhedralEdge> PolyhedralFace::edges

private

Definition at line 69 of file PolyhedralComponents.h.

Referenced by edge_sum_ff(), ff_n_core(), and PolyhedralFace().

m_area

double PolyhedralFace::m_area

private

Definition at line 70 of file PolyhedralComponents.h.

Referenced by area(), assert_Ci(), ff(), ff_2D(), ff_n(), PolyhedralFace(), and pyramidalVolume().

m_normal

kvector_t PolyhedralFace::m_normal

private

normal vector of this polygon's plane

Definition at line 71 of file PolyhedralComponents.h.

Referenced by assert_Ci(), decompose_q(), edge_sum_ff(), ff_2D(), ff_n(), ff_n_core(), normalProjectionConj(), and PolyhedralFace().

m_rperp

double PolyhedralFace::m_rperp

private

distance of this polygon's plane from the origin, along 'm_normal'

Definition at line 72 of file PolyhedralComponents.h.

Referenced by assert_Ci(), ff(), ff_n(), ff_n_core(), PolyhedralFace(), and pyramidalVolume().

m_radius_2d

double PolyhedralFace::m_radius_2d

private

radius of enclosing cylinder

Definition at line 73 of file PolyhedralComponents.h.

Referenced by ff(), ff_2D(), and PolyhedralFace().

m_radius_3d

double PolyhedralFace::m_radius_3d

private

radius of enclosing sphere

Definition at line 74 of file PolyhedralComponents.h.

Referenced by PolyhedralFace(), and radius3d().

---

The documentation for this class was generated from the following files:

• /home/www/ba/Sample/HardParticle/PolyhedralComponents.h
• /home/www/ba/Sample/HardParticle/PolyhedralComponents.cpp