CODEWITHSUNDEEP
c++cachingintegrationnumerical
Rob Simpson
Rob Simpson

Reputation: 25

Cache gauss points for numerical integration in c++

This is a question on what people think the best way to lay out my class structure for my problem. I am doing some numerical analysis and require certain "elements" to be integrated. So I have created a class called "BoundaryElement" like so

class BoundaryElement
{
/* private members */
public:
integrate(Point& pt);
};

The key function is 'integrate' which I need to evaluate for a whole variety of different points. What happens is that, depending on the point, I need to use a different number of integration points and weights, which are basically vectors of numbers. To find these, I have a class like so:

class GaussPtsWts
{
int numPts;
double* gaussPts;
double* gaussWts;

public:
GaussPtsWts(const int n);
GaussPtsWts(const GaussPtsWts& rhs);
~GaussPtsWts();
GaussPtsWts& operator=(const GaussPtsWts& rhs);
inline double gwt(const unsigned int i)
{
    return gaussWts[i];
}
inline double gpt(const unsigned int i)
{
    return gaussPts[i];
}
inline int numberGPs()
{
return numGPs;
}
};

Using this, I could theoretically create a GaussPtsWts instance for every call to the integrate function. But I know that I maybe using the same number of gauss points many times , and so I would like to cache this data. I'm not very confident on how this might be done - potentially a std::map which is a static member of the BoundaryElement class? If people could shed any light on this I would be very grateful. Thanks!

Upvotes: 0

Views: 454

Answers (1)

wxffles
wxffles

Reputation: 874

I had a similar issue once and used a map (as you suggested). What I would do is change the GaussPtsWts to contain the map:

typedef std::map<int, std::vector<std::pair<double, double>>> map_type;

Here I've taken your two arrays of the points and weights and put them into a single vector of pairs - which should apply if I remember my quadrature correctly. Feel free to make a small structure of the point and weight to make it more readable.

Then I'd create a single instance of the GaussPtsWts and store a reference to it in each BoundaryElement. Or perhaps a shared_ptr depending on how you like it. You'd also need to record how many points you are using.

When you ask for a weight, you might have something like this:

double gwt(const unsigned int numGPs, const unsigned int i)
{
    map_type::const_iterator found = themap.find(numGPs);
    if(found == themap.end())
        calculatePoints(numGPs);
    return themap[numGPs][i].first;
}

Alternatively you could mess around with templates with an integer parameter:

template <int N>
class GaussPtsWts...

Upvotes: 0

Related Questions