Some code for my calculation can be generated automatically with different
program. For now I just copy and paste the relevant output in my code which
now looks like
>
> function UnitTriangleIntegration(f::Function,NP=16)
>
>
>> xw = [0., 0., 0.]
>
>
>> if NP==16
>
> xw=[0.33333333333333 0.33333333333333 0.14431560767779
>
> 0.45929258829272 0.45929258829272 0.09509163426728
>
> 0.45929258829272 0.08141482341455 0.09509163426728
>
> 0.08141482341455 0.45929258829272 0.09509163426728
>
> 0.17056930775176 0.17056930775176 0.10321737053472
>
> 0.17056930775176 0.65886138449648 0.10321737053472
>
> 0.65886138449648 0.17056930775176 0.10321737053472
>
> 0.05054722831703 0.05054722831703 0.03245849762320
>
> 0.05054722831703 0.89890554336594 0.03245849762320
>
> 0.89890554336594 0.05054722831703 0.03245849762320
>
> 0.26311282963464 0.72849239295540 0.02723031417443
>
> 0.72849239295540 0.00839477740996 0.02723031417443
>
> 0.00839477740996 0.26311282963464 0.02723031417443
>
> 0.72849239295540 0.26311282963464 0.02723031417443
>
> 0.26311282963464 0.00839477740996 0.02723031417443
>
> 0.00839477740996 0.72849239295540 0.02723031417443]
>
> elseif NP==3
>
> xw = [0.16666666666667 0.16666666666667 0.33333333333333
>
> 0.16666666666667 0.66666666666667 0.33333333333333
>
> 0.66666666666667 0.16666666666667 0.33333333333333]
>
> end
>
>
>> s = 0
>
> for j in 1:NP
>
> xi = xw[j,1]
>
> eta = xw[j,2]
>
> w = xw[j,3]
>
>
>> s += f(xi,eta)*w/2
>
> end
>
> return s
>
> end
>
>
>> function ctrquad(f::Function,x::Array{Float64,2})
>
>
>> alpha1 = 1/(1 + norm(x[:,4]-x[:,2])/norm(x[:,4]-x[:,1]))
>
> beta1 = 1/(1 + norm(x[:,6]-x[:,3])/norm(x[:,6]-x[:,1]))
>
> gamma1 = 1/(1 + norm(x[:,5]-x[:,2])/norm(x[:,5]-x[:,3]))
>
>
>> f_(xi,eta) = begin
>
>
>> x_ = x[:,1]*(-eta - xi + 1 + eta^2/beta1 + eta*xi/beta1 -
>> eta/beta1 + eta*xi/alpha1 + xi^2/alpha1 - xi/alpha1) +
>> x[:,2]*(xi*(eta*(alpha1 - gamma1) + (alpha1 - xi)*(gamma1 - 1))/((alpha1 -
>> 1)*(gamma1 - 1))) + x[:,3]*(eta*(gamma1*(beta1 - eta) - xi*(beta1 + gamma1
>> - 1))/(gamma1*(beta1 - 1))) + x[:,4]*(xi*(eta + xi - 1)/(alpha1*(alpha1 -
>> 1))) + x[:,5]*(-eta*xi/(gamma1*(gamma1 - 1)))
>
> + x[:,6]*(eta*(eta + xi - 1)/(beta1*(beta1 - 1)))
>
>
>> e_xi = x[:,1]*(-1 + eta/beta1 + eta/alpha1 + 2*xi/alpha1 -
>> 1/alpha1) + x[:,2]*((eta*(alpha1 - gamma1) - xi*(gamma1 - 1) - (-alpha1 +
>> xi)*(gamma1 - 1))/((alpha1 - 1)*(gamma1 - 1))) + x[:,3]*(-eta*(beta1 +
>> gamma1 - 1)/(gamma1*(beta1 - 1))) + x[:,4]*((eta + 2*xi -
>> 1)/(alpha1*(alpha1 - 1))) + x[:,5]*(-eta/(gamma1*(gamma1 - 1))) +
>> x[:,6]*(eta/(beta1*(beta1 - 1)))
>
>
>> e_eta = x[:,1]*(-1 + 2*eta/beta1 + xi/beta1 - 1/beta1 +
>> xi/alpha1) + x[:,2]*(xi*(alpha1 - gamma1)/((alpha1 - 1)*(gamma1 - 1))) +
>> x[:,3]*(-(eta*gamma1 - gamma1*(beta1 - eta) + xi*(beta1 + gamma1 -
>> 1))/(gamma1*(beta1 - 1))) + x[:,4]*(xi/(alpha1*(alpha1 - 1))) +
>> x[:,5]*(-xi/(gamma1*(gamma1 - 1))) + x[:,6]*((2*eta + xi - 1)/(beta1*(beta1
>> - 1)))
>
>
>> f(x_,e_xi,e_eta)*norm(cross(e_xi,e_eta))
>
> end
>
>
>> UnitTriangleIntegration(f_,16)
>
> end
>
>
>>
>> ### Here is some testing
>
> x1 = [0.,0.,0.]
>
> x2 = [1.,0.,0.]
>
> x3 = [0.,1.,0.]
>
>
>> y = Array(Float64,3,6)
>
> y[:,1] = x1
>
> y[:,2] = x2
>
> y[:,3] = x3
>
> y[:,4] = (x1+x2)/2
>
> y[:,5] = (x2+x3)/2
>
> y[:,6] = (x1+x3)/2
>
>
>> area = ctrquad((x,exi,eeta) -> 1, y)
>
> I would like to put the code from ` f_(xi,eta) = begin ...... end` outside
of this file as it is automaticaly generated. Trying to do it with
`include(...)` failed for me as it checks variables. What sugestions do
have splitting automaticaly generated code from the one which is human
writed? Also are there anything I can do to improve performance?
>
>>
>
