rafavzqz / geopdes Goto Github PK

Dear Rafa,
Thank you for your previous help.
I would like to implement some strain computation of Cantilever Plate,but I cannot find it.Like the following image.

I only find the ex_kirchhoff_rectangular_plate.m in the folder geopdes-master\geopdes-master\geopdes\inst\examples\elasticity.
The plate is the rectangle plate.If the GeoPDEs source code has the above function,please tell me about that.
Would you give me some guidance,please?
Thank you so much.
Deeply grateful,
Jim

Dear Rafa
Two different posteriori error indicators have been introduced in your paper entitled "Algorithms for the implementation of adaptive isogeometric methods using hierarchical splines" for the scalar field. Are these indicators applicable to the vector field approximations like elasticity problems?
I mean without the exact solution at hand, can we estimate the residual error for the vector field using the above-mentioned error estimators?

Hi Rafa,

Is there anything so far developed with PHT(polynomial spline over hierarchical t-meshes) spline based isogeometric analysis in GeoPDEs? If so then please kindly let me know. Thanks in advance and have a nice day.

Best regards,
Jim

Dear Rafa
I appreciate the time you spend guiding me.
I am writing to ask if you could help me with how I can translate control points between different hierarchical levels?
suppose I have control points of level one, how can I calculate the control points of (T)HB space?
I have written the following code to calculate all control points of each level, is there any way to directly calculate the control points of each element separately?

%%
coeff_lev = cell(size(hmsh.nel_per_level));
CntrlPts = cell(size(hspace.space_of_level));
% Control points of first level
CntrlPts{1} = geometry.nurbs.coefs;
% Control points of other levels using subdivision matrix
for i = 1:numel(hspace.Proj)
CntrlPts{i+1} = CntrlPts{i}*hspace.Proj{i}';
end
for ilev = 1:numel(hspace.space_of_level)
[~,jj]=find(hspace.Csub{ilev}(hspace.active{ilev},:));
coeff_lev{ilev} = CntrlPts{ilev}(:,hspace.active{ilev})*hspace.Csub{ilev}(hspace.active{ilev},jj);
end

Best Regards
Hamid

Dear Rafa,
Thank you for your previous help.
I would like to implement Scordelis-Lo roof problem related computation.I have read the documents in the appendix and GeoPDEs code libraries this week.However,my code level is very limited.
I feel a little embrassed ... Would you please give me a simple example like how to
compute Scordelis-Lo roof problem in isometric analysis with the GeoPDEs source code?If the GeoPDEs source code has the above function,please tell me about that.(Sry,I can't understand this code very well...)
Would you give me some guidance,please?

Deeply grateful,
Jim

Hi
I solved Laplace equation on a unit square using GeoPDE, the source term is chosen such that the exact solution is given by u(x, y) = atan(25(x - y)).
After ten steps of refinement, the final mesh is given as follows:

It should be noted that the solution has a singularity in the diagonal of the domain but the calculated mesh does not show explicitly the singularity.
Would you mind explaining what the problem is?

Best Regards
Hamid

Dear Rafa
Thank you for your previous help.
I would like to implement HB splines related computation.I have read the documents in the appendix and GeoPDEs code libraries this week.However,my code level is very limited.
I feel a little embrassed ... Would you please give me a example like how to
compute and plot the HB basis functions with the GeoPDEs source code?If the GeoPDEs source code has the above function,please tell me about that.(Sry,I can't understand this code very well...)
Would you give me some guidance,please?

Deeply grateful,
Jim

Dear Rafa
I am writing to ask if you could help me with how to get the elemental connectivity matrix of a hierarchical mesh?
Best Regards
didadidadada

Dear Rafa
Thank you for your previous help.
As you mentioned the estimator for linear elasticity should read
h_Q (\int_Q (f + div(\sigma(u)) : (f + div(\sigma(u)))^{1/2}
I think displacement field is the error estimator in above equation.
I would like to know how we can make one of the components of stress or von Mises stress as our error estimator?

Best
Hamid

I am studying about the order of convergence in Isogeometric Method for the simple elliptic problem. We know that the size $h_K$ of an element $K$ in the physical mesh is related to the element size $h_Q$ of the element $Q$ in the parametric domain. My question is how can we find $h_K$ in GeoPDE ?

Dear Rafa
Thank you for your previous help.
I would like to know if it is possible to plot the THB basis functions with the GeoPDEs source code?
Of course, the hierarchical mesh can be drawn using hmsh_plot_cells code, is there any way to draw the hierarchical space?

Best Regards
Hamid

Dear Rafa,
Thank you for your previous help.
I would like to implement some computation of stiffness matrix,but I cannot find it.
If the GeoPDEs source code has the above function,please tell me about that.
Would you give me some guidance,please?

Deeply grateful,
Jim

Hi!

First and foremost, thanks for providing this package!

I am currently working with truncated hierarchical B-splines, and I would like to assemble the mass matrix for a space corresponding to a grid of the following type:

However, I was wondering what the best way of constructing such a space using your package might be?
I've tried using the build_hspace_from_cells function, but this requires me to label each element individually, and for deep nesting this may be infeasable.

Do you happen to have any convenience-functionality lying around that I have missed, or is this something that is fairly easy to do in Octave/Matlab? My Octave/Matlab experience is limited.

Thanks in advance!

Dear Rafa
I encounter a problem when I solved a linear elasticity problem with your GEOPDE package. The computed displacement of the load point is totally different from the result obtained from other distributed IGA package, e.g., igafem-nguyen. The problem setting and the corresponding boundary condition are as follows:

problem_data.geo_name = nrb4surf([0 0], [3 0], [0 1], [3 1]);

% Type of boundary conditions
problem_data.nmnn_sides = [];
problem_data.press_sides = [];
problem_data.drchlt_sides = [];
problem_data.symm_sides = [];

% Physical parameters
E = 1; nu = .3;
problem_data.lambda_lame = @(x, y) ((nuE)/((1+nu)(1-2nu)) * ones (size (x)));
problem_data.mu_lame = @(x, y) (E/(2(1+nu)) * ones (size (x)));

% Source and boundary terms
problem_data.f = @(x, y, ind) zeros (2, size (x, 1), size (x, 2));
problem_data.h = @(x, y, ind) zeros (2, size (x, 1), size (x, 2));

% Exact solution (optional)
% uxex = @(x,y) sin(2pix).(sin(2piy));
% uyex = @(x,y) sin(2pix).(sin(2piy));
% problem_data.uex = @(x, y) cat(1, ...
% reshape (uxex (x,y), [1, size(x)]), ...
% reshape (uyex (x,y), [1, size(x)]));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree = [2 2]; % Degree of the bsplines
method_data.regularity = [1 1]; % Regularity of the splines
method_data.nsub = [64 32]; % Number of subdivisions
method_data.nquad = [3 3]; % Points for the Gaussian quadrature rule

The calculation process are changed into:
fixeddofs = [sp.boundary(1).dofs(1:2:end-1) sp.boundary(3).dofs(end)];
int_dofs=setdiff(1:sp.ndof,fixeddofs);
loaddofs=sp.boundary(1).dofs(end);
rhs(loaddofs)=-1;
u(fixeddofs) =0;
% Solve the linear system
u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs);

The final cauculated results is u(loaddofs)=-2.8568 which is tollolly different from the result obtained by the package of igafem-nguyen where u(loaddofs)= -117.9371 under the same problem configuration.
Can you help me solve this problem?
Thank you very much!

I am trying to solve the quasi-linear elliptic problem $-\nabla.(a(u))\nabla u=f\text{ in }\Omega$, $u=0\text{ in }\partial\Omega$, where $a(u)=1+u^2$. I am using Newton's method to solve this problem. In Newton's method while computing the Jacobian Matrix I need matrices $(a(u^k)\nabla u_i,\nabla v_j)$ and $(a'(u^k)\nabla u^ku_i,\nabla v_j)$, where $u^k$ is the initial guess for the iteration, $u_i$ and $v_j$ are basic functions. Here $u^k$ will be in vector form. So, how can I compute these matrices?

Hi Rrafavzqz !
I want to construct a hierarchical vector space from a hierachical scalar space obtained from your package. But, I can not get the right stiffness matrix and the matrix operation is singular if I construct the hierarchical vector space as the following code. Moreover, if the level of the hierachical vector space equals to 1, the matrix is right and the matrix operation is no longer singular. Could you help me to fix it?

1. for ilev = 1:numel(hspace.space_of_level)
2.       msh = hmsh.mesh_of_level(ilev);
3.       space_scalar = hspace.space_of_level(ilev);
4.       scalar_spaces = repmat ({space_scalar}, 1, msh.rdim);
5.       spaces = sp_vector (scalar_spaces, msh);
6.       if ilev==1
7.          sp_hier.space_of_level= spaces;
8.       else
9.           sp_hier.space_of_level(ilev)= spaces;
10.       end
11.       clear spaces
12.   end
13.   for ilev = 1:numel(hspace.space_of_level)
14.       cumsum_ndof(1) = 0;
15.       cumsum_ndof(2 : hmsh.rdim+1) = cumsum (cellfun (@(x) x.ndof, sp_hier.space_of_level(ilev).scalar_spaces));
16.       if ~isempty(hspace.active{ilev})
17.          sp_hier.active{ilev} = repmat(hspace.active{ilev},1,msh.rdim)+repmat(cumsum_ndof(1:msh.rdim),numel(hspace.active{ilev}),1);
18.       else
19.          sp_hier.active{ilev}=[];
20.       end
21.       if ~isempty(hspace.deactivated{ilev})
22.          sp_hier.deactivated{ilev} = repmat(hspace.deactivated{ilev},1,msh.rdim)+repmat(cumsum_ndof(1:msh.rdim),numel(hspace.deactivated{ilev}),1);
23.       else
24.          sp_hier.deactivated{ilev}=[];
25.       end
26.       sp_hier.active{ilev} = sp_hier.active{ilev}(:);
27.       sp_hier.deactivated{ilev} = sp_hier.deactivated{ilev}(:);
28.       Id1=zeros(size(hspace.Csub{ilev}));
29.       Id2=zeros(size(hspace.Csub{ilev}));
30.       sp_hier.Csub{ilev} =[hspace.Csub{ilev}, Id1; Id2,hspace.Csub{ilev}] ;
31.   end

Best!
didadidadada

Dear Rafa
Hi
First, thank you for the time you provide me.
I have solved Fichera corner example, and I could capture the singularity and the results for three steps of refinement are as follows

Now I want to compare the results of the energy norm with GeoPDE, would you mind telling me how I should provide the input data?
My input data is as follows

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_fichera.txt';

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides = [];
problem_data.drchlt_sides = 1;

% Physical parameters
problem_data.c_diff = @(x, y, z) ones(size(x));
problem_data.grad_c_diff = @(x, y, z) cat (1, ...
reshape (zeros(size(x)), [1, size(x)]), ...
reshape (zeros(size(x)), [1, size(x)]), ...
reshape (zeros(size(x)), [1, size(x)]));

% Source and boundary terms

radius = @(x,y,z) (x.^2 + y.^2 + z.^2);
uex = @(x,y,z) (radius(x,y,z)).^(1/4);
problem_data.f = @(x,y,z) ((-3/4)*radius(x,y,z)).^(-3/4);
problem_data.h = @(x,y,z,ind) uex(x,y,z);

% Exact solution (optional)
problem_data.uex = uex;
problem_data.graduex = @(x, y, z) cat (1, ...
reshape (x./(2*(x.^2 + y.^2 + z.^2).^(3/4)), [1, size(x)]), ...
reshape (y./(2*(x.^2 + y.^2 + z.^2).^(3/4)), [1, size(x)]), ...
reshape (z./(2*(x.^2 + y.^2 + z.^2).^(3/4)), [1, size(x)]));

% CHOICE OF THE DISCRETIZATION PARAMETERS (Coarse mesh)
clear method_data
method_data.degree = [2 2 2]; % Degree of the splines
method_data.regularity = [1 1 1]; % Regularity of the splines
method_data.nsub_coarse = [2 2 2]; % Number of subdivisions of the coarsest mesh
method_data.nsub_refine = [2 2 2]; % Number of subdivisions for each refinement
method_data.nquad = [3 3 3]; % Points for the Gaussian quadrature rule
method_data.space_type = 'standard'; % 'simplified' (only children functions) or 'standard'
method_data.truncated = 1; % 0: False, 1: True

% ADAPTIVITY PARAMETERS
clear adaptivity_data
adaptivity_data.flag = 'elements';
% adaptivity_data.flag = 'functions';
adaptivity_data.C0_est = 1.0;
adaptivity_data.mark_param = .5;
adaptivity_data.mark_strategy = 'MS';
adaptivity_data.max_level = 10;
adaptivity_data.max_ndof = 5000;
adaptivity_data.num_max_iter = 5;
adaptivity_data.max_nel = 5000;
adaptivity_data.tol = 1e-5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Best regards
Hamid

Hi! It's me again!

Thanks for your previous help.

I am wondering whether it is possible to get the number of active truncated hierarchical B-splines on each cell in a hierarchical mesh? I haven't been able to extract this information from the data structures.

Is there any clear cut way of doing this, or do you happen to have any ideas on how I might go about getting these numbers?

Best regards,
Ivar