The Flex Representation Method (FRM) I Exercises

Thought I’d start a new topic for these exercises.

On question 1, I’m confused as to what I should do for delta lambda. Is that a derivative or a very small value? Should I take the variational derivative with respect to lambda as well, or just stick with the one with respect to u? Thank you!

@Rebecca you’ll notice that when you make the substitution \lambda disappears. Thus \Pi(u,\lambda) just becomes \Pi(u), and any variation with respect to \lambda just equals 0 meaning you only have one statement with respect to \delta u to solve and you don’t need to worry about variations with respect to \lambda. You were on the right track already.

1 Like

@chris, that’s really good to know. I work with that then. Thank you!

Just a point of clarification, I wouldn’t say that the variation of \lambda is 0. It is more that \lambda is now a function of u which means that the variational derivative can be passed through instead of having to stop at \lambda and you don’t end up with the second equation.

1 Like

You can in fact take the weak form from before switching over to the penalty method and just massage it into the form that you would find doing it this way

TAs,

I’ve started looking at writing the algorithm at the end of the notes and I’m running into functions that I do not recognize or understand. Can you all shed some light on them?

IE - Index mapping: I found it referenced in the notes but it never actually says what it is or does so I’m not sure how to implement it.

EG - Element to global basis function index map: I found this in the notes mentioned but it’s only written out in the example and I don’t understand how to construct an EG for my own algorithm as it’s never explicitly defined.

dN/d\xi: In the compute basis function we have another output but I’m unsure as to what it is since previously we just output N.

compte_geom changes: now we have two new inputs instead of just the quadrature point. Is this to change our quadrature point, e, to x in terms of N and dN/d\xi?

invert_map: I understand that we are trying to invert X^F but don’t know how to carry that out either.

I know that’s a lot, but I would hope that I’m not the only one with these questions.

Thanks!

2 Likes

@klbrooks These are all great questions. I think you’ll find that once you actually need to code the math you’ve been staring at for few weeks that the rubber hits the road and things hopefully start to make more sense. I will also say that I’ve just recently gone through and written the method outlined here myself so I know that there are some deficiencies that I’ll try to remedy both here and in the notes. I’m going to go through each of these pretty fast so ask more specific questions on anything you don’t understand:

IE: This one is described briefly in https://research.coreform.com/index.php/The_Flex_Representation_Method_(FRM)_I#Helpful_index_mappings but the purpose of this function is to take a quadrature point and return an element id. This whole algorithm is written as a loop on a set of quadrature points so you need to be able to take the starting point of some quadrature point (you’ll need to give it an id somehow) and compute everything you need. So if you had 4 elements, and 3 quadrature points per element you’d need a map from i\in[1,12] to e\in[1,4]. Since you start with a quadrature point, you need to be able to get the relevant element since the basis functions and mappings can look different on each element.

EG: This is a very important mapping to understand. Taking the example above, say we have 4 elements that are each degree p^2 with C^1 between each element. If you study the spline construction section then you should be able to calculate that you end up with e_n + p functions. So globally, you have 6 unique functions in your space. This is what this basis looks like:

Each of these functions becomes a unique N_A in your basis (see eq 25-30) and thus K will be a 6X6 matrix. Now lets go to the element level. Say you’re at a quadrature point on the third element (shown from 2<x<3 in the plot). There are 3 non-zero basis functions on this element so when you compute N and dNdx each will return 3 values. To fill in K I take the inner product of each of these functions represented by the double for loop on lines 33 and 34 of the algorithm. I have an a index, and a b index, but each of these is just a,b\in[1,3] and this doesn’t tell me where to put this computation in the final K matrix. Each of these local functions corresponds to one of the 6 global functions. EG takes an element and a local function index and returns the corresponding global function so that I can correclty assemble the whole system. So on element three: EG(e=3, a=1) = y3, EG(e=3, a=2) = y4, EG(e=3,a=3)=y5.

dN/d\xi: There’s nothing magic here. I plugged the equation for bernstein polynomials into mathematica and solved for the derivative with respect to \xi. This should just give you some closed form solution you can code up. There are fancier ways to take deriatives but don’t worry about those for now.

compute_geom: I’m not sure I understand this question. You are essentially trying to compute \sum_AX_A^F\hat{N}_A^F (eq. 22). you take e, nodes, and EG as input so you can pick the correct nodes X_A^F that correspond to each basis function on that element. Then you take \hat{N}_A^F as input and dot it with the vector of nodes.

Inverse Mapping: I’m in the process of updating this section in the notes since I know what’s there isn’t sufficient for you to just pick it up and do it, hopefully this should be updated tomorrow. I’ll try to convince Dr. Scott that he should spend some lecture time talking about Newton-Raphson and if you’re still confused after all that come back and ask more questions!

Hey Chris and other TAs! I’m trying to code out my compute_mesh function. I was able to create the extractions operators (ops) and the EG index map but I’m having trouble with creating the nodes vector.

I’m not sure how many nodes I should make. In class, Dr. Scott said that we need to have a node at each endpoint (\Omega^F) and that we need one node for each basis function. I currently have it coded as such for a placeholder at least.

n_N = sum( p) + n_e
nodes = linspace(omega_F(1),omega_F(2),n_N);

Where p is my degrees vector and n_e is the length of that vector describing how many elements there are. This basically just splits \Omega^F equally between all the points.

I just don’t feel like this is right necessarily. Could you shed some light on this part of the function?

It looks like you’re not correctly computing the total number of basis functions, so your EG is likely wrong as well since you would need to know how many basis functions you have in order to return the global index A. If you assume that you have p-1 continuity between elements then n_N=n_e+p.

Otherwise your computation of the nodes looks good!

P.s when you are tyring to write math here, just use the $ symbol surrounding your statement like $your math$ and it’ll look a lot nicer! (not required, but I think it looks nice)

So does that mean I need to compute the nodes of each element individually? Each element will have p number of nodes plus the end points? I’m just struggling to make the connection. I’ll send over my EG code when I can to double check with you and see where I may have gone wrong.

@klbrooks I think I might understand where you’re getting confused here. Looking at the equation for computing geometry X(\xi)=\sum_AX_A^F\hat{N}_A^F(\xi) both the basis functions and the nodes are indexed by a global index A. That means for each basis function I will have one corresponding node. If I’m trying to compute the geometry on a specific element, I’ll compute \hat{N}_a^F (note this basis function has the local index a) and then EG(e,a)=A will tell me which global function and therefore which global node should be multiplied here. So if I’m on an element, I have p+1 functions to evaluatie, but I don’t have p+1 unique nodes on each element. I have p+1 nodes that impact the current element. Interior nodes will always impact the geometry of more than one element.

The concept of geometry is really hard to understand in 1d, since we’re mapping a 1d parametric coordinate \xi to the 1d line x. Say that my geometry is actually 2d so each X_A^F has two components (x,y). X(\xi) then maps \xi to the x,y plane and it’s easier to see how each node impacts the parametric curve. Here’s an example using the same basis I used in the above post:


Notice that because I had 6 global basis functions I have 6 nodes. This is the same as what your nodes are doing in 1d but only on the x axis. As long as your first node starts at 0, your last node ends at L, and each node in between is increasing, then this line on x will only increase.

So in summary:
Compute the total number of global basis functions using the formula I gave you, and then just order that many nodes linearly from 0 to L in a vector called nodes. Then when you’re on an element and you need p+1 nodes (one for each function evaluation on that element), you’ll use EG(e,a) = A, X_A^F = nodes(A)

Ok that’s helping. Here’s where I’m still getting tripped up.

For example p = [1, 1, 2, 1]^T so there are 4 elements with different degrees respectively. Now since p is not uniform, do I need to just go with the highest degree (in this case, 2) and make the number of nodes n_e+(p+1)=7 or will it instead be n_e+p=6? Or am I making this more complicated and should be using a uniform p? That would certainly simplify things.

Thanks for all your help with this! It is very helpful! I’m just stuck on one or two things I guess.

P.S. Here is my code for computing EG:
%% Compute EG - element to global function index map
for e=1:n_e
count = e;
for a=1:p+1
EG(e,a) = count;
count = count+1;
end
end

(Looking at it now, I see that I’m assuming a uniform p I think.)

You should assume uniform p. We have only given you the extraction operators for uniform p and k=p-1. Non uniform degree is a more complicated basis to compute.

1 Like

So in trying to determine the dN/d\xi for the basis functions, I’m struggling to determine the derivatives of the basis functions themselves. Looking at the bernstein basis, I currently have it coded as such:

B(1) = nchoosek(n,0)*((1-xi)/2)^(n-0)*((xi+1)/2)^0;
for i=1:n
B(i+1,1) = nchoosek(n,i)*((1-\xi)/2)^(n-i)*((\xi+1)/2)^i;
end

Looking at the class notes, I’m not understanding how to code in the derivative to this function because I don’t fully understand the nchoosek (n;i). Then looking at equation (15) in The Problem of Approximation II notes I don’t understand the n index for B_i,n-1 (Sorry I can’t figure out how to do more than character for subscripts and superscripts).

I guess I was thinking in class we might go over some kind of function to compute the basis derivatives rather than hard code it in for each basis.

@klbrooks
Curly brackets are how you group things together in Latex for example \xi^{n-i} is produced by “\xi^{n-i}”. Also use “{n \choose k}” for {n \choose k}. As far as how n choose k affects a derivative, it doesn’t. It is just a constant.

When the notes say B_{i} they really should say B_{i,n} which means the ith basis function of Bernstein basis of order n. We say the Bernstein basis, but really it is a collection of bases, a basis for each order. So B_{i,n-1} means get the Bernstein basis of one order lower and grab the ith basis function in that basis.

The “\frac” command can also be used to nicely format fractions. For example “\frac{\xi+1}{2}” produces \frac{\xi+1}{2}

Equation 13 on the class notes where you were talking before should be
B_{i,n}(\xi) = {n \choose i}(1-\xi)^{n-i}\xi^i
Notice that as n changes the basis changes